{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "9. DiffusionReaction(PI-DeepONet).ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "machine_shape": "hm",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU",
    "gpuClass": "standard"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/jdtoscano94/Learning-PINNs-in-Python-Physics-Informed-Machine-Learning/blob/main/9_DiffusionReaction(PI_DeepONet).ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7JQXTWz35c40",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "687a5ea4-75cb-47ae-aaba-7392fe192103"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/jax/experimental/optimizers.py:30: FutureWarning: jax.experimental.optimizers is deprecated, import jax.example_libraries.optimizers instead\n",
            "  FutureWarning)\n"
          ]
        }
      ],
      "source": [
        "import jax\n",
        "import jax.numpy as np\n",
        "from jax import random, grad, vmap, jit, hessian, lax\n",
        "from jax.experimental import optimizers\n",
        "from jax.nn import relu\n",
        "from jax.config import config\n",
        "from jax.numpy import index_exp as index\n",
        "from jax.flatten_util import ravel_pytree\n",
        "import itertools\n",
        "from functools import partial\n",
        "from torch.utils import data\n",
        "from tqdm import trange, tqdm\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "from scipy.interpolate import griddata\n",
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Auxiliary Functions"
      ],
      "metadata": {
        "id": "fDhK-Kr3GyGa"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Use double precision to generate data (due to GP sampling)\n",
        "#https://github.com/PredictiveIntelligenceLab/Physics-informed-DeepONets/blob/main/Diffusion-reaction/PI_DeepONet_DR.ipynb\n",
        "def RBF(x1, x2, params):\n",
        "    output_scale, lengthscales = params\n",
        "    diffs = np.expand_dims(x1 / lengthscales, 1) - \\\n",
        "            np.expand_dims(x2 / lengthscales, 0)\n",
        "    r2 = np.sum(diffs**2, axis=2)\n",
        "    return output_scale * np.exp(-0.5 * r2)"
      ],
      "metadata": {
        "id": "UyghW-A9h4Kz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# A diffusion-reaction numerical solver\n",
        "#https://github.com/PredictiveIntelligenceLab/Physics-informed-DeepONets/blob/main/Diffusion-reaction/PI_DeepONet_DR.ipynb\n",
        "def solve_ADR(key, Nx, Nt, P, length_scale):\n",
        "    \"\"\"Solve 1D\n",
        "    u_t = (k(x) u_x)_x - v(x) u_x + g(u) + f(x)\n",
        "    with zero initial and boundary conditions.\n",
        "    \"\"\"\n",
        "    xmin, xmax = 0, 1\n",
        "    tmin, tmax = 0, 1\n",
        "    k = lambda x: 0.01*np.ones_like(x)\n",
        "    v = lambda x: np.zeros_like(x)\n",
        "    g = lambda u: 0.01*u ** 2\n",
        "    dg = lambda u: 0.02 * u\n",
        "    u0 = lambda x: np.zeros_like(x)\n",
        "\n",
        "    # Generate subkeys\n",
        "    subkeys = random.split(key, 2)\n",
        "\n",
        "    # Generate a GP sample\n",
        "    N = 512\n",
        "    gp_params = (1.0, length_scale)\n",
        "    jitter = 1e-10\n",
        "    X = np.linspace(xmin, xmax, N)[:,None]\n",
        "    K = RBF(X, X, gp_params)\n",
        "    L = np.linalg.cholesky(K + jitter*np.eye(N))\n",
        "    gp_sample = np.dot(L, random.normal(subkeys[0], (N,)))\n",
        "    # Create a callable interpolation function  \n",
        "    f_fn = lambda x: np.interp(x, X.flatten(), gp_sample)\n",
        "\n",
        "    # Create grid\n",
        "    x = np.linspace(xmin, xmax, Nx)\n",
        "    t = np.linspace(tmin, tmax, Nt)\n",
        "    h = x[1] - x[0]\n",
        "    dt = t[1] - t[0]\n",
        "    h2 = h ** 2\n",
        "\n",
        "    # Compute coefficients and forcing\n",
        "    k = k(x)\n",
        "    v = v(x)\n",
        "    f = f_fn(x)\n",
        "\n",
        "    # Compute finite difference operators\n",
        "    D1 = np.eye(Nx, k=1) - np.eye(Nx, k=-1)\n",
        "    D2 = -2 * np.eye(Nx) + np.eye(Nx, k=-1) + np.eye(Nx, k=1)\n",
        "    D3 = np.eye(Nx - 2)\n",
        "    M = -np.diag(D1 @ k) @ D1 - 4 * np.diag(k) @ D2\n",
        "    m_bond = 8 * h2 / dt * D3 + M[1:-1, 1:-1]\n",
        "    v_bond = 2 * h * np.diag(v[1:-1]) @ D1[1:-1, 1:-1] + 2 * h * np.diag(\n",
        "        v[2:] - v[: Nx - 2]\n",
        "    )\n",
        "    mv_bond = m_bond + v_bond\n",
        "    c = 8 * h2 / dt * D3 - M[1:-1, 1:-1] - v_bond\n",
        "\n",
        "    # Initialize solution and apply initial condition\n",
        "    u = np.zeros((Nx, Nt))\n",
        "    u = u.at[:,0].set(u0(x))\n",
        "    #u = index_update(u, index[:,0], u0(x))\n",
        "    # Time-stepping update\n",
        "    def body_fn(i, u):\n",
        "        gi = g(u[1:-1, i])\n",
        "        dgi = dg(u[1:-1, i])\n",
        "        h2dgi = np.diag(4 * h2 * dgi)\n",
        "        A = mv_bond - h2dgi\n",
        "        b1 = 8 * h2 * (0.5 * f[1:-1] + 0.5 * f[1:-1] + gi)\n",
        "        b2 = (c - h2dgi) @ u[1:-1, i].T\n",
        "        u = u.at[1:-1, i + 1].set(np.linalg.solve(A, b1 + b2))\n",
        "        #u = index_update(u, index[1:-1, i + 1], np.linalg.solve(A, b1 + b2))\n",
        "        return u\n",
        "    # Run loop\n",
        "    UU = lax.fori_loop(0, Nt-1, body_fn, u)\n",
        "\n",
        "    # Input sensor locations and measurements\n",
        "    xx = np.linspace(xmin, xmax, m)\n",
        "    u = f_fn(xx)\n",
        "    # Output sensor locations and measurements\n",
        "    idx = random.randint(subkeys[1], (P,2), 0, max(Nx,Nt))\n",
        "    y = np.concatenate([x[idx[:,0]][:,None], t[idx[:,1]][:,None]], axis = 1)\n",
        "    s = UU[idx[:,0], idx[:,1]]\n",
        "    # x, t: sampled points on grid\n",
        "    return (x, t, UU), (u, y, s)"
      ],
      "metadata": {
        "id": "HF71T_Koh1G9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Plots\n",
        "def plot(X,T,f):\n",
        "  fig = plt.figure(figsize=(7,5))\n",
        "  plt.pcolor(X,T,f, cmap='rainbow')\n",
        "  plt.colorbar()\n",
        "  plt.xlabel('$x$')\n",
        "  plt.ylabel('$t$')\n",
        "  plt.title('$f(x,t)$')\n",
        "  plt.tight_layout()"
      ],
      "metadata": {
        "id": "It-HA0ea2w-f"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# **Introduction**\n",
        "\n",
        "## **Functions and Operators:**\n",
        "\n",
        "### **Function:**\n",
        "\n",
        " Maps between vector spaces:\n",
        "\n",
        "*Example:*\n",
        "\n",
        "Let $f_1(x)=sin(x)$; for  $x\\in\\mathbf{R}$\n",
        "\n",
        "$$z=f_1(x)=sin(x)\\in[0:1]$$\n",
        "\n",
        "In other words $f_1$ maps $\\mathbf{R}→[0,1]$\n",
        "\n",
        "### **Operator:**\n",
        "\n",
        " Maps between infinite-dimensional function spaces:\n",
        "$$G(f_1(x))=f_2(x)$$\n",
        "\n",
        "*Example:*\n",
        "\n",
        "Derivative Operator →$\\frac{d}{d x}$\n",
        "\n",
        "It transforms a funcion $f_1$ into a function $f_2$:\n",
        "\n",
        "Let $f_1(x)=sin(x)$\n",
        "\n",
        "Then when we apply our operator:\n",
        "\n",
        "$$f_2=\\frac{df_1(x)}{d x}=\\frac{d}{d x}sin(x)=cos(x)$$\n",
        "\n",
        "\n",
        "### **Parametric PDEs and Operators**\n",
        "\n",
        "Parametric PDEs $\\rightarrow$ Some parameters (e.g., shape, IC/BC, coefficients, etc.) of a given PDE system are allowed to change.\n",
        "\n",
        "Let  $\\mathcal{N}$ be a nonlinear differential operator. Let's consider parametric PDEs of the form:\n",
        "\n",
        "$$\\mathcal{N}(u,s)=0$$\n",
        "\n",
        "where, $u$ is the input function, and $s$ is the unknown PDE's solution (also a function.)\n",
        "\n",
        "Our PDE solution operator would be:\n",
        "\n",
        "$$G(u)=s$$\n",
        "\n",
        "**Note 1:**In other words, we can express the general solution of our PDE as an operator $G$\n",
        "\n",
        "**Note 2:** Remember $s$ is itself a function, so if we evaluate it at any point $y$, the answer would be a real number:\n",
        "\n",
        " $$G(u)(y)=s(y)\\in \\mathbf{R}$$\n",
        "\n",
        " ### **Universal Approximation Theorem for Operator**\n",
        "\n",
        " $\\forall \\epsilon >0$, there are positive integers $n,p,m$, constants $c_i^k,W_{bij}^k,b_{bij}^k,W_{tk},b_{tk}$ such that:\n",
        "\n",
        "$$\\left|G(u)(y)-\\sum_{k=1}^{p}\\sum_{i=1}^{n}c_i^k\\sigma\\left(\\sum_{j=1}^{m}W_{bij}^{k}u(x_j)+b_{bi}^k\\right).\\sigma(W_{tk}.y+b_{tk})\\right|<\\epsilon $$\n",
        "\n",
        " ### **Neural Network**\n",
        "\n",
        "A Neural Network is a function that takes the form: (See:https://book.sciml.ai/notes/03/)\n",
        "\n",
        "$$NN(X)=W_n\\sigma_{n-1}(W_{n-1}\\sigma_{n-2}(...(W_2\\sigma_1(W_1X+b_1)+b_2)+..)+b_{n-1})+b_n$$\n",
        "\n",
        "So we can use 2 NNs to implement the Universal Approximation Theorem for Operator i.e.\n",
        "\n",
        "Branch:\n",
        "\n",
        "$$NN_b(u(\\textbf{x}))=b(u(\\textbf{x}))=\\textbf{c}.\\sigma\\left(W_{b}u(\\textbf{x})+\\textbf{b}_{b}\\right)$$\n",
        "\n",
        "Trunk:\n",
        "\n",
        "\n",
        "$$NN_t(\\textbf{y})=t(\\textbf{y})=\\sigma(W_{t}.\\textbf{y}+\\textbf{b}_{t})$$\n",
        "\n",
        "### **DeepOnet**\n",
        "\n",
        "Learn the solution operators of parametric PDEs → We will try to approximate $G$  (the solution of our PDE operator) by two neural networks:\n",
        "\n",
        "\n",
        "$$G_\\theta(u)(y)=\\sum_{k=1}^q\\underset{Branch}{\\underbrace{b_k\\left(u(x_1),u(x_2),...,u(x_m)\\right)}}.\\underset{Trunk}{\\underbrace{t_k(\\textbf{y})}}$$\n",
        "\n",
        "\n",
        "We want to obtain G, so our goal would be:\n",
        "\n",
        " $$G_\\theta(u)(y)\\approx G(u)(y)$$\n",
        "\n",
        "So we will enforce that condition into a loss function:\n",
        "\n",
        "$$\\mathcal{L}_{Operator}(\\theta)=\\frac{1}{NP}\\sum_{i=1}^N\\sum_{j=1}^P\\left|G_{\\theta}(u^{(i)})y_j^{(i)}-G(u^{(i)})y_j^{(i)}\\right|^2$$\n",
        "\n",
        "\n",
        "$$\\mathcal{L}_{Operator}(\\theta)=\\frac{1}{NP}\\sum_{i=1}^N\\sum_{j=1}^P\\left|\\sum_{k=1}^q{b_k\\left(u(x_1),u(x_2),...,u(x_m)\\right)}.t_k(y_j^{(i)})-G(u^{(i)})y_j^{(i)}\\right|^2$$\n",
        "\n",
        "where:\n",
        "\n",
        "$m:$ Number of points at which we evaluated our input functions.\n",
        "\n",
        "$N:$ Number of input functions.\n",
        "\n",
        "$P:$ Number of points at which we evaluate the output function → output sensors.\n"
      ],
      "metadata": {
        "id": "PhPIO4gwxTcb"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Physics-informed DeepONets\n",
        "\n",
        "\n",
        "Similar to a PINN, a Physics-informed DeepONets output functions consistent with physical constraints by minimizing the residual of an underlying governing law (i.e., nonlinear differential operator).\n",
        "\n",
        "$$\\mathcal{L}_{Physics}(\\theta)=\\frac{1}{NQm}\\sum_{i=1}^{N}\\sum_{j=1}^{Q}\\sum_{k=1}^{m}\\left|\\mathcal{N}(u^{(i)}(x_k),G_{\\theta}(u^{(i)})(y_j^{(i)})\\right|^2$$\n",
        "\n",
        "where $\\mathcal{N}$ is a nonlinear differential operator, and $\\{y_j\\}_{i=1}^{Q}$ are the collocation points ( we use them to enforce the physical constraint).\n",
        "\n",
        "So the total loss would be:\n",
        "\n",
        "$$\\mathcal{L}(\\theta)=\\mathcal{L}_{Operator}(\\theta)+\\mathcal{L}_{Physics}(\\theta)$$\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "uRXrweOQ7JRL"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Problem Setup\n",
        "\n",
        "**Diffusion-reaction system**\n",
        "\n",
        "Implicit operator described by a nonlinear PDE with a source term $u(x)$:\n",
        "\n",
        "$$\\frac{\\partial s}{\\partial t}=D\\frac{\\partial^2 s}{\\partial x^2}+ks^2+u(x)$$\n",
        "\n",
        "$$(x,t)\\in (0,1,]\\times(0,1]$$\n",
        "\n",
        "where, $D=0.01$ is the diffusion coefficient and $k=0.01$ is the reaction rate.\n",
        "\n",
        "**Note:** We will not use any paired input-output data. We only know that  the initial and boundary conditions  are zero.\n",
        "\n",
        "\n",
        "**Training**\n",
        "\n",
        "We will map the source terms $u(x)$ to the PDE solution $s(x,t)$. Hence, we will approximate the implicit solution operator ($G$) with a PI-DeepONet ($G_{\\theta}$).\n",
        "\n",
        "From our PDE, we know that for a given input function $u^{(i)}$:\n",
        "\n",
        "$$u^{(i)}=\\frac{\\partial s^{(i)}}{\\partial t}-D\\frac{\\partial^2 s^{(i)}}{\\partial x^2}-k[s^{(i)}]^2 $$\n",
        "\n",
        "So since ideally $G_{\\theta}(u^{(i)})(x,t)\\approx G(u^{(i)})(x,t)= s^{(i)}(x,t)$ :\n",
        "$$u^{(i)}\\approx \\frac{\\partial G_{\\theta}(u^{(i)})(x,t)}{\\partial t}-D\\frac{\\partial^2 G_{\\theta}(u^{(i)})(x,t)}{\\partial x^2}-k[G_{\\theta}(u^{(i)})(x,t)]^2 $$\n",
        "\n",
        "Lets call:\n",
        "$$R_{\\theta}^{(i)}(x,t)=\\frac{\\partial G_{\\theta}(u^{(i)})(x,t)}{\\partial t}-D\\frac{\\partial^2 G_{\\theta}(u^{(i)})(x,t)}{\\partial x^2}-k[G_{\\theta}(u^{(i)})(x,t)]^2 $$\n",
        "So our $\\mathcal{L}_{Physics}$ will be:\n",
        "\n",
        "\n",
        "$$\\mathcal{L}_{Physics}(\\theta)=\\frac{1}{NQ}\\sum_{i=1}^{N}\\sum_{j=1}^{Q}\\left|R_{\\theta}^{(i)}(x_{r,j}^{(i)},t_{r,j}^{(i)})-u^{(i)}(x_{r,j}^{(i)})\\right|^2$$\n",
        "\n",
        "where $(x_{r,j},t_{r,j})$ are the \"collocation points\" (i.e., where we evaluate our PDE).\n",
        "\n",
        "On the other hand we will enforce our zero initial and boundary conditions with the $\\mathcal{L}_{Operator}(\\theta)$:\n",
        "\n",
        "\n",
        "$$\\mathcal{L}_{Operator}(\\theta)=\\frac{1}{NQ}\\sum_{i=1}^{N}\\sum_{j=1}^{P}\\left|G_{\\theta}(u^{(i)})(x_{u,j}^{(i)},t_{u,j}^{(i)})- G(u^{(i)})(x_{u,j}^{(i)},t_{u,j}^{(i)}))\\right|^2$$\n",
        "\n",
        "where, $(x_{u,j},t_{u,j})$ are points form our initial and boundary conditions. Consequently, since we are working with zero initial and boundary conditions $G(u^{(i)})(x_{u,j}^{(i)},t_{u,j}^{(i)})=0$ and:\n",
        "\n",
        "\n",
        "$$\\mathcal{L}_{Operator}(\\theta)=\\frac{1}{NQ}\\sum_{i=1}^{N}\\sum_{j=1}^{P}\\left|G_{\\theta}(u^{(i)})(x_{u,j}^{(i)},t_{u,j}^{(i)})))\\right|^2$$\n",
        "\n",
        "Finally, the total loss will be: \n",
        "$$\\mathcal{L}(\\theta)=\\mathcal{L}_{operator}(\\theta)+\\mathcal{L}_{Physics}(\\theta)$$\n",
        "\n"
      ],
      "metadata": {
        "id": "J-qnbwQ5DTuw"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**In summary:**\n",
        "\n",
        "*   Let $y=(x,t)$.\n",
        "*   Also remember that $u$ is only a function of $x$.\n",
        "\n",
        "To train a PI-DeepOnet:\n",
        "\n",
        "1.   Select $N$ source terms (i.e, input functions) →$\\{u^{(i)}(\\textbf{x})\\}_{i=1}^{N}={[u^{(1)}(\\textbf{x}),u^{(2)}(\\textbf{x}),...,u^{(N)}(\\textbf{x})]}$.\n",
        "2.   Evaluate our $N$ functions at $m$ points (i.e., input sensors) →$\\begin{bmatrix}u^{(1)}(x_1)&u^{(1)}(x_2)&...&u^{(1)}(x_m)\\\\\n",
        "u^{(2)}(x_1)&u^{(2)}(x_2)&...&u^{(2)}(x_m)\\\\\n",
        "\\vdots&\\vdots&\\ddots&\\vdots\\\\\n",
        "u^{(N)}(x_1)&u^{(N)}(x_2)&...&u^{(N)}(x_m)\\end{bmatrix}$\n",
        "3.   Send the $m$ outputs of our $N$ functions to our **branch network** → $b_k\\begin{pmatrix}u^{(1)}(x_1)&u^{(1)}(x_2)&...&u^{(1)}(x_m)\\\\\n",
        "u^{(2)}(x_1)&u^{(2)}(x_2)&...&u^{(2)}(x_m)\\\\\n",
        "\\vdots&\\vdots&\\ddots&\\vdots\\\\\n",
        "u^{(N)}(x_1)&u^{(N)}(x_2)&...&u^{(N)}(x_m)\\end{pmatrix}$\n",
        "4.   Select $P$ points (i.e., output sensors) from our boundary and initial conditions → $\\textbf{y}_u=y_{u1},y_{u2},...,y_{uP}$\n",
        "5.  Send our output sensors to our **trunk network**→$t_k(y_{u1},y_{u2},...,y_{uP})$\n",
        "6. Approximate our operator by computing the dot product between the output of our **branch network** and the output of our **trunk network**→ $G_\\theta(u)(\\textbf{y})=\\sum_{k=1}^q\\underset{Branch}{\\underbrace{b_k\\left(u(x_1),u(x_2),...,u(x_m)\\right)}}.\\underset{Trunk}{\\underbrace{t_k(\\textbf{y})}}$\n",
        "7. Ideally $G_\\theta(u)(\\textbf{y}_u)\\approx G(u)(\\textbf{y}_u)=0$, so compute the error → $\\mathcal{L}_{Operator}(\\theta)=\\frac{1}{NP}\\sum_{i=1}^N\\sum_{j=1}^P\\left|G_{\\theta}(u^{(i)})y_{uj}^{(i)}\\right|^2$\n",
        "8. Select $Q$ points (i.e., collocation points) inside our domain → $\\textbf{y}_r=y_{r1},y_{r2},...,y_{rQ}$\n",
        "9. Compute the error related to the underlying governing laws (i.e., PDE)→$\\mathcal{L}_{Physics}(\\theta)=\\frac{1}{NQ}\\sum_{i=1}^{N}\\sum_{j=1}^{Q}\\left|R_{\\theta}^{(i)}(y_{r,j}^{(i)})-u^{(i)}(x_{r,j}^{(i)})\\right|^2$\n",
        "\n",
        "10. Compute the total loss→$\\mathcal{L}(\\theta)=\\mathcal{L}_{operator}(\\theta)+\\mathcal{L}_{Physics}(\\theta)$\n",
        "11. Update our NN parameters (i.e., branch and trunk) to minimize  $\\mathcal{L}(\\theta)$.\n",
        "12. Repeat the process until $G_\\theta(u)(x,t)\\approx G(u)(x,t)$."
      ],
      "metadata": {
        "id": "RyGFoStnSGGs"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Data Generation"
      ],
      "metadata": {
        "id": "FvY9PrxfylwH"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We will randomly sample 5000 different functions $u$ from a zero-mean Gaussian process with an exponential quadratic kernel with a length scale: $l=0.2$."
      ],
      "metadata": {
        "id": "phC4C1GNDrvA"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# GRF length scale\n",
        "length_scale = 0.2\n",
        "N = 5000 # number of input samples"
      ],
      "metadata": {
        "id": "R-lQbo2kYcvA"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Generate a random function"
      ],
      "metadata": {
        "id": "gupNZ3GKDSDu"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "Let: $$(x,t)\\in (0,1]\\times(0,1]$$\n",
        "Lets create a random function:\n",
        "$$UU=f(x,t)$$\n",
        "We select $P=300$ points from our BC and IC (i.e., $100/side$)\n",
        "$$y_u=\\{(x_{u1},t_{u2}),(x_{u1},t_{u2}),...,(x_{uQ},t_{uQ})\\}$$\n",
        "We select $Q=100$ collocation points to evaluate our PDE.\n",
        "$$y_r=\\{(x_{r1},t_{r2}),(x_{r1},t_{r2}),...,(x_{rQ},t_{rQ})\\}$$\n"
      ],
      "metadata": {
        "id": "QRK8Q3jNtewv"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Resolution of the solution (Grid of 100x100)\n",
        "Nx = 100\n",
        "Nt = 100\n",
        "# Select the number of sensors\n",
        "m = Nx   # number of input sensors\n",
        "P_train = 300 # number of output sensors, 100 for each side \n",
        "Q_train = 100  # number of collocation points for each input sample"
      ],
      "metadata": {
        "id": "k1COdMm5uID4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "config.update(\"jax_enable_x64\", True)\n",
        "# Numerical solution\n",
        "key = random.PRNGKey(0)\n",
        "(x, t, UU), (u, y, s) = solve_ADR(key, Nx , Nt, P_train, length_scale)"
      ],
      "metadata": {
        "id": "zle1RAx7dFoN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "x = np.linspace(0, 1, Nx)\n",
        "t = np.linspace(0, 1, Nt)\n",
        "XX, TT = np.meshgrid(x, t)\n",
        "plot(XX,TT,UU)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 369
        },
        "id": "pxfo2SFadNyK",
        "outputId": "b9dacbd9-960e-4a42-f1f6-0fd99c365355"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 504x360 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdwAAAFgCAYAAAARq8j7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2de7AlV3Wff+ucc6+EA+bhwUKWRFAK4SADNjAREGIjG8kRVBmVXzwcYmGDlZjg8gO7UEop7II8hCnjR0wZjzHFww9sEx7jQo4iZDAJhYiG8DAaGWuQMQwIhGQgEMHc81j545xhVq99e/XjdPfp0/f3Vd2a23fv3nt333Nnn/2d1WuLqoIQQggh7TLa9AAIIYSQgwAnXEIIIaQDOOESQgghHcAJlxBCCOkATriEEEJIB3DCJYQQQjqAEy4hhBDSAZxwCSGEkA7ghEtIA4jIg0XkRhH5ooi8TkT+i4j8XMN9/G8R+Y4m2ySEdIcw0xQh6yMirwJwtqq+UEQeDODDAB6uql9bs92TAH5AVT8kIs8E8CxV/eEGhkwI6RiucAlphssA/Nnq++cBuL6ByfYQgHMAHF/96CiA7xWRh6zTLiFkM3DCJWQNRGRXRL4M4NEA/lxE/hrA0wD8lanzqyLydnP8ShG5SUR2g3YfDuDTWP6N3iMi9wCYAfgggH/ZztUQQtpksukBELLNqOqeiDwJwLtV9RwAEJEvAPi4qfYKAHeIyGMBPAHAFQD+haruBe2eEJFfBPA9qvqs0z8XkdsAfGcLl0IIaRlOuISsz3cB+Ig5fgCAr5w+UNV7ROTXAbwBwP2xnGy/XKLd78Tys2DLVwCcu95wCSGbgEqZkPXxE+4XAdzP1fkQltr536vqp2u2i1W7X6ozSELIZuGES8j6fCeyE+NHATzi9IGIPBrA72C5wv3JMg2KyAjAo5CucB+JdBImhGwBnHAJWR8/4V4P4CkAICLnAfhzAP8WwAsBPFpELrUni8jrReT1rs37rL5Gpt7ZAB4P4MZmh08I6QJOuISsweoRnQcC+Bvz4zcCeLqI3B/LyfdVqnpUVe8F8EoA/8k1cwGA99kfqOr/A/AaAMdXz+ICwA8AeI+qfrb5KyGEtA0TXxDSAiLynwHcpaq/UVBvF8vV8WNUdVpQ9wMAnq+qH2tupISQruCESwghhHRA60p5lVf2LhHZ9125LPktETkhIh8Vkce1PSZCCCGka7r4DPf1WD7on8fTAFy0+roay2hOQgghZFC0PuGq6nsB/ENQ5UoAb9QlNwN4gIjwwX5CCCGDog+Zps7DMmfsaU6ufnanrygiV2O5CsY/Ejz+n561Gn6Vz6GH8pE1P3snhPSAD04Xd6vqg5tu9woRvbvmuR8EblDVyKxuhD5MuKVR1SMAjgDA4fvs6LELH7AsmFeYfGaLep1X6aML5jWvg3RD314vhLSE3PmVv2+j3bsBHKt5rgCHmhxLU/Rhwv0Mls8hnub81c8IIYQcZMY1P/Xs6YKkD4kvjgL48VW08hMBfFlVE51MCCHkACEAxlLvq6e0vsIVkT8GcCmAQ6uMOb8MYAcAVPU1WGbieTqAEwDuBfATa3caaeNI9dV9V9SUPuziXVlP3/mVgpr24LHNr1eyJlJ/hYt5oyNpitYnXFV9TkG5Avh3bY+DEELIFiEAJv1drdahD0qZEEIIGTx9CJraHEW6qqzCrKK9mlJkbenVTSu82QHQxpu+x4RsA4I1lHI/OdgTLiGEkP7S4wCoOnDCJYQQ0j9knaCpfjLMCXdifklRxLL/ZbalhsvWXUcTt6Epu9C7fderjIwuT99/l2S7oFImhBBCOoJKmRBCCGkZKuWe4t8FNaUBbbu+TftCKFJpZeu2dR1VaEMjd5LQY4vUL9UrIeVoaYUrIlcA+E0AYwCvVdXrXPlZAN4I4PEA7gHwLFX95Lr9DuvtAyGEkGEgWMbj1PmKmhUZA3g1lnuxXwzgOSJysav2fABfVNWHA/h1AK9o4pI44RJCCDlIXALghKreoap7AN6M5b7slisBvGH1/VsAPFVE1l5uD0MpRyTvdozOizQxkFV/ke6tEu1cqW4FxRx91lFFYUap1Orq5nWiwUv3EYy7b7q5qd8VIUNG1tqI4JCI2N39jqy2dwX234P9Ce78b9RR1ZmIfBnAt2C5a2Bthj/hEkII2U7qB03draqHmxxKExzsCbdoBRkFO4XvvCqsYCq9oGo+z1ulj2iFVTeRuF8ZtxF5WMUUtEFjgXpb9CkPV+OkTdp7DrfMHuyn65wUkQmA+2MZPLUWW/TXTQgh5MAgbo/b5vbDvQXARSJyoYjsAng2lvuyW44CuGr1/Y8A+MvVznZrcbBXuIQQQvpLCyvc1WeyLwJwA5aPBb1OVW8VkZcBOKaqRwH8PoA3icgJAP+A5aS8NsOccKPnZ6O0j41p4prtFLZb9tw1VN94bPpvKtioC/W4YVkzdscHQbeO/UU7+hasRs7Q5Mcsey1t9i5o7eMgVb0ewPXuZy81338dwI823e8wJ1xCCCFbzvAyTQ3ragghhJCecrBXuAUZSTLKuSlNXEWvemVXVtE19dxrU5o0Uo+NRfdWqNuF7p1UGZChix2auqLmLSDbxvYp5U1xsCdcQggh/YTb8xFCCCEdwRXulrFOyr9IOZfWzZ513rE1oHTXUbhtpCTchAouiq4tSxtRuNSwZNv42rSddrk9HyGEENIRnHAJIYSQlmHQ1MBYZ8P3oghnS2397Gng3V5bCRqqaNq6Krapd7uNXXMzzfSOg5C0wzOwldQgoFImhBBCOmJgK9xhvX0ghBBCegpXuJa676bWiXauRBsbt3eQlCLps+Z5fYtS3gRd5CceisYb2OqoVfr4O+dzuIQQQkgHSKmt9rYKTriEEEL6CVe4JKGzd2Edv/iqvDr8Vodt0DcVvInt53p2C1pnYCucRmjsI6qew8eCCCGEkC7gY0GEEEJI+3CFOwC6eMfUVuKAPr/46ureTWjZpjh4fz3t0OfXdZP0ebXWx7ExSpkQQgjpguEp5WFdDSGEENJThrnC3fS7ok33v00ctKhbkrLNfy+b1uHbfO+KEACjYX3cMMwJlxBCyPYzsDcUnHAJIYT0D0Ypk04Z2Ls70gAD+w8IQP9e55u+x13cj01fYymGFzTFCZcQQkj/4AqXEEII6YgRV7jEMzDt8Q0G9u5ybYbye+7b77UP97WLe7KtqnhTv58B7hbUg1c6IYQQMny4wiWEENI/mNqxp3Sianrwi98mvdKH+1WWPt/XTd/HTdwbqtd22m3qmrt8Tfb5b7MGw5hwCSGEDAsRBk0RQgghncAVLmmNoURKFo7hACQWSPrs+Jq3VcsC7Yx9qFq2b+NpkgF+htvJ1YjIFSLycRE5ISLX7FP+UBF5t4h8SEQ+KiJP72JchBBCesxI6n31lNZXuCIyBvBqAJcDOAngFhE5qqrHTbX/AOBPVfV3RORiANcDeFjbYyuklXfaW7QqSPoYyAp8m6+jz8E+Q121bXo8dftf5/o3baGA1XO4XOFW5RIAJ1T1DlXdA/BmAFe6Ogrgm1ff3x/AZzsYFyGEkD6zgRWuiDxIRG4UkdtX/z5wnzrfJSLvF5FbV1b2WaUuZ62RleM8AJ82xydXP7P8CoDnishJLFe3P7NfQyJytYgcE5FjX5gt2hgrIYSQg801AG5S1YsA3LQ69twL4MdV9TsAXAHgN0TkAUUN9yVo6jkAXq+qvyYiTwLwJhF5lKpmZlVVPQLgCAAcvs+ObmCc9di0kgrb7LHi7vPYkjZ79jvug07d9GvgICjcKn1u4rrWYXNBU1cCuHT1/RsAvAfAS2wFVf1b8/1nReQuAA8G8KWo4S4m3M8AuMAcn7/6meX5WL5LgKq+X0TOBnAIwF0djI8QQkjvWEsPHxKRY+b4yGrBVoZzVPXO1fefA3BOVFlELgGwC+ATRQ13MeHeAuAiEbkQy4n22QB+zNX5FICnAni9iDwSwNkAvtDB2AghhPSR9Va4d6vq4dymRd4F4CH7FF1rD1RVRSTXporIuQDeBOAqb2T3o/UJV1VnIvIiADcAGAN4nareKiIvA3BMVY8CeDGA3xORn8cygOp5qtq9Mu6Dhsu0w/HE7XSsDBvrf8P3se/qs25/fb+usn2sdR0lz62yctxkpHBLj/io6mV5ZSLyeRE5V1XvXE2o+5pWEflmAO8EcK2q3lym304+w1XV67EMhrI/e6n5/jiAJ3cxFkIIIVvA5h4LOgrgKgDXrf59h68gIrsA3gbgjar6lrIND+shJ0IIIcNhM4kvrgNwuYjcDuCy1TFE5LCIvHZV55kAvgfA80Tkw6uv7ypquC9RygTox4PqjOxs/jxgmCp2m/too90q/beldNvQ1ptKgrGhKGVVvQfLmCL/82MAXrD6/g8A/EHVtjnhEkII6SH9TtNYByplQgghpAO4wm0CPpzfXD+dRItusd7soo/SCrUH96PPuncT96fs/rFdRX+vg4D74RJCCCGd0IdNFBqEEy4hhJD+IcIV7qDYpqjTKn32rY++aeKDqFf7fM3brHurtBNNHl2Mp8rkFbXTpWIeWNDUwZ5wCSGE9BMBlTIhhBDSCVTKpDJtKLou+ti4pt2AAuuzwi1qt6nXQFO/n7I6cJ372NQ119W9XUQX1+2/Kf2+Ka0rgsXAlPKw3j4QQgghPYUrXEIIIb1DASyolA8oW6Vie6wzq9TtSic20U5RG01FfzeilNdQwXX1b5/V6yb6rHuP1zm3jfNaZmhKmRMuIYSQ3qEimG9yL94W4IRLCCGkl3CFu+1s01ZX2xQF20kfHdy7PvdfVN61Ct5E5G8b6reoz7L3tWhy6Fr3+vO2TSkLoPwMlxBCCGmXZdAUV7iEEEJIuwzwOVxOuHXpW1KGLiJ2Nx2h2pqm7VrFVtB326yCNx3524XC7ZvurXtf/XnheLqZBIf4WNCwroYQQgjpKVzhEkII6SVUykOmKVUy1CjYztVrBf29EaXcMzXcNxXcp8jfKskkutDE/j62oXuD31XRRGZV7qYmPRXBXIYlYTnhEkII6SVc4RJCCCEdwAl329i0wl3n3K4jZH27m77mpvpoKplEF/q7impsSzdvq6Ztqm6VexXVTT5WKPd7b0r3+rSItm5RQglbd3NKmYkvCCGEkA7gc7iEEEJI+8jwnsPlhBvRhQruQht3oa3bqhudV/YeVIl0bSuatymNPqo51rp6ta0kEGV/P01F/hb9x5255vJ9ZNVreU0b6V6P16plda8fT6b/gujfqI+uNK8CWMiwVrjDevtACCGE9JSDt8Lt22ovOq+pNItt1G0r+KrKNdcNaKobbNTF6nedQKS6z7a2sTItqls2MKnoGdSagUjRatTXtavRKqu9tM/81XAY/CT5deuOJ+lf8seattPdOo2f4RJCCCEtoyL8DJcQQgjpgvnAPsMd5oRbZXeasJ0WnhetpCU7UMFdBEa1oYmL+mhD6VbRzU1p2r61U1ax+7qRCq6kZYNAoJYCkaLnXqPgo0gpx/3Hujfu48y5frKyZVpFIQdKu02GuFvQMCdcQgghW44kbwy2HU64hBBC+ocwaGpYVNoYfQN6ddPaeB2lXfc5yyrXUTbStUo7dZ9tLYoYrpsusa42XmenmpoquPYzoRVSEFaJ/I3Oqxv5G+rmNSJ9rf71ZZH+jXSvH888Gmuwkkx+Bx2tOpfP4Q5LKQ/ragghhJA1EJEHiciNInL76t8HBnW/WUROishvl2mbEy4hhJBeshhJra81uQbATap6EYCbVsd5vBzAe8s2PHylXGlHoEDRrRPNm2lnHU3bgO5trI8K6rPKeOruarNOBHETqQyLFHJZbVx0X8vq7yhhRME9L5ugoSgqOFaxJfsoiJAtG4lcGO2cUbrtRP5aPZqWmWsuSkphjiPdG52XqOhEP+fX7SyQSWRTqR2vBHDp6vs3AHgPgJf4SiLyeADnAPjvAA6XaXj4Ey4hhJCtQ5G++ajAIRE5Zo6PqOqRkueeo6p3rr7/HJaTagYRGQH4NQDPBXBZ2UFxwiWEENJL1ljh3q2quatOEXkXgIfsU3StPVBVFRHdp94LAVyvqielwhgP3oRbJWFF3XaiiN3G1HQLZf64bnRxUZ9ldW+Rbo70dxeRv3WVcpVdbVrIHRzlCvblVRJEzCbjbLuRQm0oCURUt2+Rv1qyblqWfw/mgbZOrhE1+6gQ0dwkbe4WpKq5q1IR+byInKuqd4rIuQDu2qfakwB8t4i8EMB9AeyKyFdVNfq89wBOuIQQQvqPSGdbATqOArgKwHWrf9/hK6jqvzr9vYg8D8DhoskWYJQyIYSQnrJYBU5V/VqT6wBcLiK3Y/n57HUAICKHReS16zTMFa6lbhKGohzItdtpQyl3FF1c9rrqamJfvukI4o6SSVht25QKrhslXLQ1XFk1nKroblWwL49U8Cw8r5nIX6t+92snE+3syqZWWyPoI9DNSR8F42mLTW1Ar6r3AHjqPj8/BuAF+/z89QBeX6ZtrnAJIYSQDuAKlxBCSC/Z0HO4rcEJ19KUpm0tgrgBxe115u44v+46ySSa0MY7wdh83bailO0YNpBMYu50axPbyCVlYR/1EkQs+wwiiMdWi7ajgufjcW5ZFU1rVfB85O5VC5G/cxnnlqV9joIydx2wffixBe14peyO20JFmEu5DiJyhYh8XEROiMi+kVwi8kwROS4it4rIH3UxLkIIIf1lQ0FTrdH6CldExgBeDeByACcB3CIiR1X1uKlzEYB/D+DJqvpFEfnWtsdFCCGkvywzTfV38qxDF0r5EgAnVPUOABCRN2OZq/K4qfNTAF6tql8EAFXd70Hj8tRVr030V7XPKlHKVv9W0caRJvZKuamkFLbdpqKLI8W8TjsNRBAnkbY1I4h9O1XUcEZNJ+0ESRcCNTwdZ9vRSCkHati3UzYq2JelUcLrq+B0PPkqeOr7N5IwiaAOVPAs0Mbzgqhgq5wj3evLrMaewf1eK7TTlVIGlXItzgPwaXN8cvUzyyMAPEJE3iciN4vIFfs1JCJXi8gxETn2hdmipeESQgjpAypS66uv9CVoagLgIix3aDgfwHtF5NGq+iVbaZV8+ggAHL7Pzn75LQkhhAwARYer6Y7oYsL9DIALzPH5q59ZTgL4gKpOAfydiPwtlhPwLaV6aEwN14z8bSpKOYkYrhAZbTVplf4jTRtp40RF11S6VaKLExVcMhJ6UnBfg+jesokmokjf5XF+H1a3VkkmMfOatmRe4XlwHpBVqjMfNV0hd3DZKGGvDadBBLNfvcwkv24UFeyVbtmkEFVUcBQJHGnaaRAxvDw+Uz5Hft1EKdsoZXfeXPPH43XzTN3fU4v0OQCqDl0o5VsAXCQiF4rILoBnY5mr0vJ2rPYfFJFDWCrmOzoYGyGEENIJra9wVXUmIi8CcAOAMYDXqeqtIvIyAMdU9eiq7PtF5DiAOYBfWqXXIoQQciAZXtBUJ5/hqur1AK53P3up+V4B/MLqq12qbM9XVhuvE6W8W1IFA7GKtedWiTyuUresCq5St0rkcZXo4p18Tet16zSom8mjO/F6t5wm9uVeE1ttXRRBnNXP+bo10s1RxLA/TvRzoF5nLhI4yh08D6KLbd0kKjjRvUavJhHE+f3P/HiihBElo4LbUsGzoK5XwbZdn0t5YerONB6rbXdeULctNpVLuU36EjRFCCGEnEHSz963HU64hBBCegdXuNtIkTYuW7dI95atG0Ui77pfRxT5G7VTJfI4qhtFBVfJcxzVjSKIfeTxju/D6t58hVolP3HUTl1NDGT/46iyHV0SGV0zgjiKGI6SSUyD/qOI4aSdKAeyj5jNbDEXRxeXjRJOlHKkaf11mSQRScQugrE2pIKnTumq2nsnuXXV9WE1su/fX9d0YcbjHsD07baHJPdi2xn+hEsIIWQr6XMSizpwwiWEENI7qJSHRlHCjLLaOEqKkdStkCCiSkRzFBUcJcVoSylntrUL+oy0cZGmtbrXjadsdDEQJ56w+jlJihGcl+pno/NqamIgq38TjR3mIDYJPPw1BnmGfeRx2bzCgIuaDnRzFDFctI1cpHSnmUQTvo983Tv112HLgojductPvEi0cX7kb5WoYKt7q6jg+SL/Xs0Wrs9FoL9dXVKegz3hEkII6S1M7UgIIYS0jDLxxZZQJYI4Oq90LmWvgoPkElWSWfhI5EgN7waaNtoqz6vhsmq6Sjs1tXGiaV0foe61OZC9hnSKWYNoZ6t4IxWdROzWTDzhE02kEcT5urfsdnSRJgayqngeRAUnEcOR0q2QTMJq4ihiGPDRva5/BOq1phqOEkZE0cT+3ETTRmWLCu3Mo3asbnb3Y+6ilk35dO6uubMoZa5wCSGEkNZRYdDUsFjr2doqQVOjcmWFuwVJft1oJ5+6K9NohVu0A48t38lf7fkVpV3Fps+25q+Go1VstBk7kF3FNhXsVHcVG+2c48uTazYr17qrViC7co1WrVUCmqK0h8mK0q7Gkd+mr5usMG07boebuitVf549TlaNi/xVo28nEzTlVpTTRTCehV/FmnbceXalqm478dk8Gk+84m4T/2z1tnOwJ1xCCCG9ZPkZLidcQgghpHX8Z+/bzvAn3CLdW/bcurrZn5s8vxo8o+vrlg2MCgOYKijlaHee5Dyn6DIBTU5tmXO97rUqOApS8nWj1IqFujfQxrMgJaJtZ51gpzD4KtDG0abqyfO7JTXxst18xV1WEwNZVRypYF9mlW7Upq/rg5SmmeCi8s+2hs/I1gxg8n0kAVWmXa930wCncmo4bce0MXPX6LSxfQ7X153Nu5sEhxalPKyrIYQQQnrK8Fe4hBBCtg4FHwsaFtFzt/sd551bpK1HgTaONHHSTslI5CopGasoZat/3XUsEjV85rju87Ne/ab6OdgBJ2qnBW1c+GyrbadmdPHyXKNpK6RLzEYXx8/PRtrYKt2i3XEy7VR47jVSwdHzq+tEEJdNe+ifSa3ybKut69W0VbpeBddVwwvXh63rtbCvO59ZpRzXbQ/hhEsIIYR0ASdcQgghpGUUaWDgtsMJtyxh4ouCBBqRfi6riX25jza2x1EkchWl7JNbmLJo5x6gKBI5X/dGCSt8usRs4gunxEzd5LwgajlOSuFU4yj/OiJtHG3G7nVvEqVsypuKLo7SJXpNnNnxpmCjdKuKo2QSyXklNTGQ1biRJi5Ku2hVsU8YUTaCONHN8/w+1OteU9dHN89m/rpM3VlwHYEK9pp4NnXHplzcWHem3U2CXOESQgghLaOQZFvFbYcTLiGEkF5yYBNfiMjDVPWTLY6lHbrIgRxFJSd91izz7VaJLq6ilIMcyGUjj5d1jRb0+ZJLbvgeKWTfThR5nGhrH/lr9XOyqXs5bbzODjyzTHSxu8ZAG/u6cVIK039BfuIoujijm9fIT2y1sc8VnC3L18S+D7+JellNDGRVcZRoItLESZmLIM7kQPYqOIggTpTyzF4XsmVTc+/ceWqOvRaeuPtxtimf+Lru+G60x9CUcpX1+lv9D0TkiQ2OhRBCCNkoIvIgEblRRG5f/fvAnHoPFZH/ISK3ichxEXlYUduFE66IPFNErgNwPxF5pEjmrfORshdBCCGElOV04os6X2tyDYCbVPUiADetjvfjjQBeqaqPBHAJgLuKGi6jlN8H4GwALwDwKgDfLiJfAvBZAF8rcX73RAkrWumvgpqOklsUJb6oq7ijPipEIpfNgezPjRJYRBu+pwkr8qOW60YeL4/tWJ3eNGVRXuMoQYU/N9LGSXRzC9o4yl3sz420caSJfV2vm6dBEoi6+YnT7egkKMuPKI7K0oQV5hoDTQxkVXGkiRMVnKjqM8fq6lpVfJZTvzum7niaHdt4njnMaOMDqJSvBHDp6vs3AHgPgJfYCiJyMYCJqt4IAKr61TINF064qvoZAG8UkU+o6vtWnX0LgIcB+JtSwyeEEEIqoJB19sM9JCLHzPERVS1rZM9R1TtX338OwDn71HkEgC+JyFsBXAjgXQCuUdX5PnW/QemgqdOT7er7ewDcU/ZcQgghpCprRCnfraqH8wpF5F0AHrJP0bWZ/lVVRHSfehMA3w3gsQA+BeBPADwPwO9HgxrGY0FNKeR1IojLtlulj6aSWwSRx0lOZNOujyC2x5FCBpzSDRJfVFLBQQKLSD/vJUo5GGsFbVx2O7zlcb4KzuRH9ro3SFJRVxsn59XUxkVJKcpq46aiiyMVXJSfuGziiUQp25zDvv9EMZt2nNItq4mBrCrecXWtKvbnjczaq0gTh0p5b/sTX6jqZXllIvJ5ETlXVe8UkXOx/2ezJwF8WFXvWJ3zdgBPRMGEO6ynigkhhAwCxfIz/Dpfa3IUwFWr768C8I596twC4AEi8uDV8fcBOF7UMCdcQgghvWRDUcrXAbhcRG4HcNnqGCJyWEReCwCrz2p/EcBNIvLXAATA7xU1PAylXJd1VHQm8reCbo5yKVeJdo624IvG4zVx7W318lWwr1s3Etkr5FM72ZerjTausnVelQQWpzK6ufuEFX7rOlt3L9rWLtDGifqtqY2ThBUtaONU/dbTxslWeYE2biq6OMpPPHFKeXcWKNwqurdu2Z6vm182CsOCtp9VjNJT9/n5MSyf1jl9fCOAx1Rp+2BPuIQQQnqJQg5uakdCCCGkS7h5wbbjNW3ZutF5RSrYJ5vI1A1UcBSJXDe5hdPEXg3bnMRVttULI5prRiJHCtmfWyUH8nTk9XO5BBbRVnmnRl5Tl9fGZSOPgaxGrrId3rxCwopMdLEr27OJHrxubkEb+/zEdbXxqb2CCOJgW7umoovvUzO6eOwTWEyDMqN/q5VlDjPlO143u7pt4vNabzsHb8IlhBDSexRYJ/FFL+GESwghpIdIZoelITD8CbeuQm60bpT4okKUchSJPIr6KJfMAvCRyOUTXyT6ObMFX6CCg2QWfou7upHIkUJelhud6FTw3ujMn4hXytkcyPkKGQBO2SQZQeTx1J3nFfNeJtq4mYQVewt3P4zuTXIQ20joIK+xP3fPRf6W1cZRggo/hiiX8Z5TwW1o40TT1tTGVhnv265RvLtfz441ii62xxPXR1h3L67bFqc3LxgSw59wCSGEbB8WCCIAACAASURBVB+abl6x7XDCJYQQ0ju4wh0adbVwUVlUN4pYLopSzvQRREL7LfdMWZKwIolalvy6RtOqK/N1y0Y0N7etXvlkFlEkslXIQPkcyEkSiiCBRbXI4/IJLE6ZsiqRx1OXsGJqFPM6eY73Zua6am6H55NQTBPda/VzfqKJuW8nSFIxdxHNZbWx3/KurjY+y2viCkkpymrj9Lyobn4ZqcbBnnAJIYT0FgZNEUIIIa3TSF7kXsEJNyJUwxW0cSYSeQ1tXDZ/c4Ut91I1bHReoIJ9O14Nz02EsYqvm7+tnY1M9gp54aOEa26rF+VE9n/gZSORu8iBDGSjlMPkFk732kjkKOcxkNXIp2ZO1ZuyJGLY9zm1unezCSv2Trn+Z17Tnik/KynLz0G8a9qNElT4cp9MomwZAJz1tfW18XpKGZ2gYOILQgghpBMYpUwIIYR0ADcvGDJVttUrW1al3SjRxX7HmbL8MSwy29jlK2RfN9LPPio5OvbJLSJtbCOTfeKLqRu7VdVJXTtWr2mDnMhRJLKPLs7o5kAhA9mEFlEkcqSQ/blJnmOjjaOt85LEF04bZyKavSa229EF0cVAVhsniS/MsT8vU+b07qkkgYUZT4XI490ogjhIWHHW18vr5nB7PKd0d0y7UV7jZXm5CGJftvu1/P4jbVxUty2olAkhhJAuUKFS3grWWY020n/NlXK0gvXn+hVtJu1j/iq2KLVj2bpRkBTgVp8+faQZXxK0lVn9+pWxq4uorlndFAVfZQKj8vuIdvlJgq2CFI1NbQ7vV7HZlWl+YFS0ol2eG6xMK+zkk1kNu3bsuaf2nB0oueE7kA2G8itcm4bRPxO764Ko6m7cvnOqfLBTtFtP2UAofxyvYvPbKXq2NrvCdX3c280kuFzhdtJVZ3Qw+wAicoWIfFxETojINUG9HxYRFZHDXYyLEEII6YrWV7giMgbwagCXAzgJ4BYROaqqx129+wH4WQAfaHtMhBBC+g8TX1TnEgAnVPUOABCRNwO4EsBxV+/lAF4B4Jc6GFM5ulDTVZ7DzZTln+cVrsUrXI9Xqtl2y5dl1LTTzYsg2CmrgvM1NZDVyF5p22Anr3t90FR2PF4b5z+jmwl2crLIR1daNZwEcQW7/Pj/cKIUjVa/Rc/WRgoZyD4XGwVGRQoZKB8YVSUlo3+etuxOPpFC9sfps7b5KjjUxD7AKnpG1mhkq4X3rVsyoCltp9x5Sd1T7n58DZ0wxKCpLpTyeQA+bY5Prn72DUTkcQAuUNV3Rg2JyNUickxEjn1htmh+pIQQQnrDYpVtqupXX9l40JSIjAC8CsDziuqq6hEARwDg8H12BvZxOiGEkNMomPiiDp8BcIE5Pn/1s9PcD8CjALxHlmrvIQCOisgzVPVYqyOLoombbKdA455pp6BeqJ/zy+yzt9Fzt748iiD253nda9VslL4xVcqjfb8HUqVbVk37duauHat/vRq2qtqrYDse32aihm20s1fBkH2/B1I1nGnHRSLPM7rZjcfo30RTBykavdKLNof3kdHRrj/2+dooEtlHHleJRB7XfEZ291R+O1EkcqSQ/bGPRK72jGwbUcqunXvL120NFX6GW4NbAFwkIhdiOdE+G8CPnS5U1S8DOHT6WETeA+AXW59sCSGE9JrFghNuJVR1JiIvAnADgDGA16nqrSLyMgDHVPVo22MghBCyXVAp10RVrwdwvfvZS3PqXtrFmCpTabP6NZJblK1bM2I4rVuvzEcFhwk0fFIMm3YxUMGppvaRyOV0b6Sil8f5kchZbRwp5QL9jXyla/9T8ZHHiRpe5LeTGY/TvVYN+43ik/EsbDv5Y/VJCZI+bTvBLj/qxmO18cyp34WLlbQpGqPN4cX1ESnlkbvmKBLZlok7r8rG8bbPakkp6pa5Pk5VqQtSk40HTRFCCCEJOrzHgjjhEkII6SUMmiLNUkVVR9HOJSOhI01cVJ7Nsxz357Vt2TKf7CJ7Xn6fXj9H/fmkFBkV63Vz8Ki6j0wO+7AqNujfv6NP29m/TcCr6Xxt7KOJfTuZsXr9bI595LHXvRk17ftY2O9dH/P8/r0aHpnjsS8z7fgI5rSdM9973Zzpw2njbARzpgji7kdWG+crXK+0R3O4uvn62Z6bnme+DxSyPzdqp00UwhUuIYQQ0gVD27yAEy4hhJDeoZramG2HE66lUiRyhRdC2bpVIpgDirSxpW7u5KRu1E6oe4PzClK0lU2gEanfon40iHa2ZYkyTbbnK6d7k7FFZX48mQjiWBtnxpZoY1OWtLN/f0Cqf+1xpKbniYq2bWTLJk63ilWfrq7Vxl6Ljt1xpGIl0KvjTJnklvlzfTtly9ZqZ1avDx+V7Ou2CT/DJYQQQjqAn+ESQgghLaOgUt4+msqX3DcqaOOiLfnKUkUxh+3UVczBeUndQBPX1dh1y4Dyaqwos06on4P/nCIVnNYtr58jNNjQK9LNliia2FOkdKN2fERxth2jmytsUlYUbZxXN9W7+fcn6iPur9oxaYYutucjhBBCqqHLx4LqfK2DiDxIRG4UkdtX/z4wp96visitInKbiPyWSPGKgBMuIYSQ3qFYWpI6X2tyDYCbVPUiADetjjOIyD8H8GQAj8Fyt7t/BuApRQ0PXyl3QRXV2lROZkeVyOSIpvSz35KvDkUKOUpKsU1EOjpSzE0ldm8qMKWpz9uq7BDjk13k4dV0G0RaukkOku7dUNDUlQAuXX3/BgDvAfASV0cBnA1gF4AA2AHw+aKGOeESQgjpH7qx7fnOUdU7V99/DsA5voKqvl9E3g3gTiwn3N9W1duKGuaESwghpHesuT3fIRGxe6ofUdUjpw9E5F0AHrLPeddmxqCqIpLkuxKRhwN4JIDzVz+6UUS+W1X/ZzQoTrikFlHuYtINY/f/wFBM42h05rqKkpbMTd1JlMBk5P/PbP71q2t8wrEY2/EVfJQyzjuv+Nz92ugzfuvGCtytqodz21W9LK9MRD4vIueq6p0ici6Au/ap9oMAblbVr67O+QsATwIQTrjD+BCMEELIoFAsH2er87UmRwFctfr+KgDv2KfOpwA8RUQmIrKDZcBUoVLmhEsIIaR/qGC+qPe1JtcBuFxEbgdw2eoYInJYRF67qvMWAJ8A8NcAPgLgI6r650UNUymTWoi2v43HqIM+uiLzMZC7rJH/gSHVxvX+M4mC2EcV+hhb3evzCjttm4mc94kVTN00oP3MeV4F+2h8a0a9XrX/8U6m2R4Wrs8oqH5u2vWJJuZmAP4/Uz+ehUnCEyldXxa1E52btGMGuHC/j7TPcmUAhvNZxgpVvQfAU/f5+TEAL1h9Pwfwb6q2zQmXEEJI71BsLEq5NTjhEkII6SXcLYik+D3EEPgivxfZToW6k/y6o4Z2apbMtdQPZRybdmbO1zWlikc408fchSNEmnbkU9GY3Mr+PGvLwjZdmcK/Js70kWhicxiVAVk17e9jRtNq/nX4pxy8Cp6b8khp+3YiLTtK1LA9cpp4lFcvq3cBYGI1rbvltk+vYf2flv37me76snwVbM/zZf6zxFEQXTzbzc+lHCndmRur/bzC6+8outn270nG4xVyW0pZ9/mvdcvhhEsIIaR3UCkTQgghXaDcno/0hQoKOauJs45u7PTqLHhSbGTbGWc9V6JpTTuR3vTn2eO5lFfRSZn5Ox073btwx1YHj107M6tw3XmCQL264dgxqG8n6MNHEI9MwhGV/HaS8dgy9ytWd832tvv7OjYqcu5eg/7ztp3Jmd+lOjdso413Ju73Y9r12+HN3Nhn5lzxkdDmJTrdif9eMlvw+T6Cc6229So6CUfPkL+t3my3KJlFfpKMbDtB9wn5fRbmbt6r0k95FMIVLiGEENIFDez80yuY+IIQQgjpgOGvcL3b2/QYdoJ6DUUaRxHLVaKZRz5E0LhIXzYfr6F/bRdmfHMX2Zrq5zN9JprY1B15JecOrcZNlK6JNlYfTRto4rGLUl4YTyvJWPPbmYhvx0TMuuvYMXXVJ4jQ/I8VvBreGefXtbp3J8nj68ZqNPJ4EtXN9rFT4TVqxz6K/raCqHGg/pZ3XmO7Undc9rqKzpNSZWnCjKgPT762ThJffKWgqbpoc1tQ9oXhT7iEEEK2DkYpDx3/cF74PG1DK2fbZ/RMLuAeSvMP5+WX2VWtf0PuV7z2WIOy5NlNXzd4XnRhxjoKntH16SP96tMGOC3cSnCidiXoVlBupWxXriO3C5Lt06+iJ5nngN2K0gUtLfTMnfcbzu8EO5f7umNb1907G/w08dFGmbG4/oOVqo8SnYyjD9XyV8PxCtOtjCd59fbjzLl7blW/G56b7fOUqVt3Jx115/nH5m158Eh9Qt20i2lZ/srYP6ecaccHtQWpJZuGz+ESQgghbaNrbc/XSzjhEkII6R1UygeZdVIy2rpRQIg/L6qbuJb8gCarmIuCpkL9bJSl7yPRv0bbjhf5gT9JekIb7OTVr9O9E5Njbu7KrP6d6NyV+eCnM/2E+tn97VsVvZMo0yxWOU/iSJsMfqyBicXYKtUkz6H53v3Vn5r5Xm3l7FjnUf7GgKnTvbaZU079nmUubJrs6uPbscfZynumnYnX717/GgXvA84yzwyfctrcPpfstaxrx5anKSvN9+6537FLu3jW12z/XgXnp7OMyiZ77uOSPfu9q5s8J9wSus+nfFsOJ1xCCCG9g4kvCCGEkC5QQOeccIdDY5HGrh3/PGAYXVxTMfs+zbHfmWVs2vGay6thqyUXM6diTVkSzbvI17+pNjbRzm4Xc6ubdwp8Ukb1BQrVpxycBNubLJx+jh55XFRQw7adPf8rr/B/yil7EG2aLt5ZBue5/wWm9rXk2rH6ede9zn2X9nUoTun6Z6wtM/N7lVl+2bLdM8c+ZeV0dqZs5iJtT+3lp5qce6Vs9Ovc7wg0Nf27Te4TNWzqznekdNnu17PtmsfPE71r9e/EpVzMauL8srSdfN1MqnGwJ1xCCCG9RMHPcAkhhJBO4Ge428A8ULhNtROWVdDEZSOY/blOe2XOdXp1lNl9xemhmYsgNuov2almka+td7waNu1M5oHCTSKPg4QNPhLZJoEI3glXUr+unUxQapBxbyr+dxf8abl25qVT/mWZ+hSR5ncyC9Jnik9SEuxeNJ8HZUGyEwCYSv5raZYpy45vahTuxKng2Uxc3TPteG08MtrYRzfPvX6emF1/nNLNaGPXh9W/YxdSbc9bltvzyuvmqYtS3tmzffoyDcoiTZwftZyWoRuUiS8IIYSQTpCaK9weZNDfF064hBBC+odm9ymuQvJoeU84eBNuGyq4KCmGbddHZ5bVxP44KTPtuBfpxEQb+yhlfxzWnZVXw14xlz0vJvgz8p8cBCrqVJVPGUw7iTYemT+fJPI4O1axyUdc5Zndrcj1cSrZkP5MOz7ZyDzKQ53ZnD6OGLaJSqaSn5d7Oo+Vss27PHMO1547cX8/E1N35vqY+oQNRvF63TwxiR523N/EzH2Usmt2tp+6dmZGze55pW008mTmxuqU7sQc1y3z5UmZ1cZTzS3zuyOto5/bQgAEKcG3koM34RJCCOk/KhgxaIoQQghpnyrxjtsAJ9yISDFHSSh8aJ3VWV43R5rYq2qreJPQTvPK9PluzXgmTgv7jePtlny+bqZJr6JLKmSP32C9Gme0rU+ucWpiXtoFXUxtggavUE3ZyHnrsblXe/48X9eo4alTuqNAN4vTtHNzMWO/kbwZw1hdmTlOlLKra8t3nG7eM7rX5yeeLbw2tnXzdbPfAnBu2pk6Fb3r1LAtn/uyWbkyAJgbxez7nBuN7HXz3Ghbr5tPueOdTJIMF8G8yNfE0XHdMj+R7ThNbMe346KS/dhJeTjhEkII6R2i2TciQ4ATLiGEkF7CoKltI8kNFkQXJyHotm6ge6OoZCAbmRzp5kgTA9lksVEEsz/PDsVpwJ2gbhVOFVcpxcjoTa+Jpz6ZwyR4+Zog4enYK1yf29noVqc+p3bbQ58IJEjsMHWqfGLqTpxuntnt6JzuTesGargB3ezbrVI2dauRXfNam/o+jbadOhWtZkvEmWszjXbO18aziY129udlfz92PItF9m8io6Z9lPI8KEv0s73mTBG+NrMq2L0GfPR1oI3tc6s7SeINW4bcMiD7OE6ipruKUlZgxM0LCCGEkPapm/iir3DCJYQQ0jtEgTGjlKsjIlcA+E0sHe1rVfU6V/4LAF6ApQj8AoCfVNW/b2UwVfIsz/OTSVRS0/Zh/T33CrKKdxxoYiBUxdk2fQINc+zaaOsFMLV5l32OW6OG/bZ+04WJ2PW5eoNc0377t5G5lyOnMxdO91rlPFGngo3unLkPlKY2X3SgkAFgZpRuXd28PDeqm6+b1ezPN010sxurjs33TsUGSnnmo5ZN+cJtkWgTWMx9dLNaTezK3McDZ5t2I/2cRB57/WzVsBurrZtENwdlPvH+zNT1STrmVls7pT13ivlec67OAhVcU0UDcUS119HtMbzncIMdNZtBls9ZvBrA0wBcDOA5InKxq/YhAIdV9TEA3gLgV9seFyGEkB6jy4xYdb76SusTLoBLAJxQ1TtUdQ/AmwFcaSuo6rtV9d7V4c0Azu9gXIQQQkhndKGUzwPwaXN8EsATgvrPB/AX+xWIyNUArgaAh04aeK+QRB5X2GPCamSnucJ2fS7lKlHCUeKLUb42rkIXLwgxGnk0yc/V67eYS/ID22jWcbadqelj7FTn3NW1ynnmfj9i6tbVzUBWOUe62W8l6KOWd41GniG/7o7TzWr2BJy5NlPFbPSz08Zn23bUKW2vnzVfPy+MUvZRyrauz53s9fPCjidIvJGc5xX3PF9x28QckYr2kdDJJ0IlE3FEKhrI6mhd+HbMWJ2K/rp5iSZ5p11d+3JOlPKsG80rGF7QVBcr3NKIyHMBHAbwyv3KVfWIqh5W1cMPbmLCJYQQ0k9WQVN1vtZBRH5URG4VkYWIHA7qXSEiHxeREyJyTZm2u5i1PgPgAnN8/upnGUTkMgDXAniGqjb1WCchhJAt5PRuQXW+1uRjAH4IwHtzx1YuNimhC4N4C4CLRORCLCfaZwP4MVtBRB4L4HcBXKGqd1XuoYoKrttupImjnMdANmp5XLMMqK+Kw3zN+ZHRExdBbCOKZ+6aZZGvhicuimFmNPLURyJPbHSxj1LOj2ge+y3ezNvcSDcDWeW868dqleHItRNEN/ucv/Zcr5utUp6PvEKW3LqJbjbHC6+tzXtrf95Zro+51daBfla3BWGS3MLU9fp5bvr012zbUad+vbYOI6EX+Zp44e9rkL/ZtptEVJvfsx9rFBmdKm1zzV4Tu9dSNkkHsnWDJB02+tor5Spq+t6usj9tKPGFqt4GABJvHfqN2KRV3dOxScejk1pf4arqDMCLANwA4DYAf6qqt4rIy0TkGatqrwRwXwB/JiIfFpGjbY+LEEJIv5FFva8O2C826byikzp5DldVrwdwvfvZS833l3UxDkIIIdvBMvFF7RXuIRE5Zo6PqOqRb7Qt8i4AD9nnvGtV9R11Oy1imJmm6qpg/8u1x1XKoqhln/iiKSJt7HM9l27TKWVzHbsFSSkWRn9Pd/J1s0+KsTDb/E1n2fMmE398ZnwLp38y2nqcr5uX5xr16X53U6Ojd5xC3TVlPumDV7pWP+86NT0LVHCkn8/212yvw20XaO+PV8pzX9ckop4lEcxGS7pIaJ9Aw2rbpE/TbpJowpQtkK+bl+cG7VjdDK+bK7SzCFRwoL8TpWz6jKKki9rJKG73f4+NcE7zR8u+9Yra8Qk93Kcn+DzaY41nau9W1dyApwYWeaVikzzDnHAJIYRsNaLoc6apwtik/eCzNYQQQsgKEflBETkJ4EkA3ikiN6x+/m0icj2QH5tU1DZXuJZQP1co28uqxyzRVn4V/Invc7ekNk40ujveKZnQwynkiY+gNhq5im6eGxWcRDfPXZSwUc4LFxk+MXV3vZaNtLGPEh7n617bznScr3CX7YxzyzIRzG5sZyX6146nfCS0VbM+KjiJ2JX8sVo1rMF5QFZHz/1Yzd+BJmX5KthHTWe0tdPNGiTpmPuEGg2oaX8/vLZeZD7NKp/Aw+vfTJ8+ijsTme0ikU3dSFMDLqK6YIvENpWytPQJXISqvg3A2/b5+WcBPN0cJ7FJRXDCJYQQ0j9U1gma6iWccAkhhPQO0X5vRFCH4U+4USRyYXRxsD2fPd5zfXq9Gypm23+Fd3N+Kz9LpI39NfqIZnvsE3hY3TxzfwkugtiW++hie50TF4lso4vnbqw2KhkAdkcmmtb1sZieKfO6d8dd1+7Equl8bezLbLtJxLBXzKbPSD/7PM+pNg4imgP9XVZFL+vavMsuJ7MpS9pJop3zo5Qz1xycF0VJ+7phO0mUdH5ktFfKmT7cR0JWxfqkHF4/23aT3NL2Xi3y1bjvcxpoYz8e23+VSOhIW7dNA1mjesXwJ1xCCCHbhwJCpbwFZFZ0QTBREgjl306Zd31JYJRZ4e2606IVb7TaLXpx2TH4Z2vt2P0K256XlPnccOa6fB9Rmc8Ybst9QJVZpY128letvg8fGBWthm3dXb/6jdrxKzFTdx60E61ol8f5AVZ21VjYjn2etuZKuXCFO7J1m1kp+7rzIDAruzKNV/yZFJHJKjoIvgpWzmkfZqx+tWl/HyhYmZYM8Cp69tiuOJMgruC5YFvXB0klq/PMajz/ueQ2Eay/EUHfGOaESwghZLvhZ7iEEEJI+wg2s3lBmxy8CTejhgtSO4ZlRsH4dI2RYo4Co4oCmqxi9Sp4HOhv22fwvGylur7Mj93W9cFXtq4LmsrUPZVtc+TGmtHP/r4ahZsEVLnxWB0daeO6Ktofp30EmnYSaFoXxDUPgqYy1+GfrQ30syYpK43uLUhnGenvTPBVoGkjTe3rehVrnzdO+yifhjJzP4J2vKZOgqYa0Na+vG7wV6q78zV2pKZJNQ7ehEsIIaT/aGc7/3QGJ1xCCCG9g0FT20iyOX0QwRzVTSKP7YHTqV4x22dmQ6XsVXCF52nL6mbff6LDx+Xq+rK6Sjkq87r7VIV2jE6c+HaSPs2x16Q2Etnr3cwzuu55yArtZHRzBTXt+7DtJhoyo8bjZ1tthHMaJZ1/zYkmtVrf9amBRm8qorqstvbnJmkxA20dPZdcJWWmVcO+LNLYyXVFajyIhE77CO5dV0p5QxvQt8nwJ1xCCCHbB6OUCSGEkPZZRilvehTNMowJN1GoZXfOCRQygKwqdmWZyOMgKtiX+7J5UOY1bUb3VtDNVv35Mt/nNNDGo0BNR+0mkdAllXJU5svrlvlyp3tH5jr8rkPZPoJ7hVjpZpRykMADyEY4R8k1kvMiFRzo8LCPNSKqMyq2IW0dqeko2tofRxp7HozV73oUaewoojuKxF62G/QZavQoaru+4m4NKmVCCCGkfYa4wh0VVyGEEELIugxzhRslt4g0cfL+Y1GuLIoKBmKlXLbM9xMp3EplwXuuSL1WUcpexZYda3SeP3fTarpAW49s/uhIP0cR1K5uJaVcM6K6rrYuaqeuto7qesoq7WSsFZKGlFXaQP3I7FA/11TcVfrweOXdGgyaIoQQQjqAEy4hhBDSPgJh0NTWMy+piX353CezCHRikqM5UsFRUowgX3GkV6vka67STqSf60ZGF40n02YHajrps4VI6CLFXbKdUaCtw6hp36anZCKQKNoaqKax8/uIFXLZpCG+HU8TirsooUhem0D5qO2kz0j9Rgk8At1dVLczuMIlhBBC2kc44RJCCCHdwAl3G4i0sdVlPrrY150HSSDK6mZPUb7kTFmgn5O6keIOxla3bnSNRXXLqukq6nWdsZZN6FFXd/t2knNrqum67URtFNS1GjtpJboHUdIQT0mlvTwup6aLEopY6kZtR/0DcfKR6Lx4PEH/FfS3p0qEd1sIE18QQggh3TC0FS4TXxBCCCEdMPwVbpE2zhDtdhzoZk8VpRydl9SNFLc9r0I0ce12GtK9Ud0i9Vl27Gso1LC/prR1XptFdYs0dl67hVHKJeu2pb+DNkdecUf9R9q0IPd17hgq/D6a0t9R3bLR3vsdR2OL6vpo8NZg0BQhhBDSPoxSJoQQQjqCE+62kyhmS6BKqkQi+xeJVUSVVHSQ+CI8NzqvSkR1BTUdjqdKuw3p1qb6iM6re3/Wib6ue14VvRr1UUVx1+2jisaOzqtSt4XxJIlJ6vbRhhovaidosyiJSFMwSpkQQgjpCK5wCSGEkLbhZ7g9JdTEAV7rVGnHvhCK9GpZjVz04ooinPPa9A0XjTVsx5Z10M462rqpdrrQvWXPW6vPhvRq2bEUnVu2nYIEDa38DppS40m7Danyts8rONdHircFg6YIIYSQjhjahMvEF4QQQkgHHOwVbl0VXdROFeUSvYNLlHfJt3uh6qsQ+Ry2WzCWJnRi3TaL2q3UTkPX0fX9GGw7FX53VaJpmxh7Ux+BrKPRq4ynbJtF42kJKmVCCCGkCxQYzTY9iGahUiaEENJLRnOp9bUOIvKjInKriCxE5HBOnQtE5N0icnxV92fLtH3wVrjzkpG+VfA6ZtOqug1NDVTTgnUjo8P+qoy1qT63tY+GdOI67Q6lj0qv+z73UfO+dpTowrNBpfwxAD8E4HeDOjMAL1bV/yMi9wPwQRG5UVWPRw0fvAmXEELIVrCJCVdVbwMACfb9VdU7Ady5+v4rInIbgPMAcMIlhBCyXay5wj0kIsfM8RFVPbL+qFJE5GEAHgvgA0V1hz/htqGQ2+qjSCs1oarbiqj2lM0fXYW6Sruw3QZCGeoqwaSdhsIq+jaepN2ejY+vgaCdDYX6rDfh3q2q+37+CgAi8i4AD9mn6FpVfUfZTkTkvgD+G4CfU9X/W1R/+BMuIYQQYlDVy9ZtQ0R2sJxs/1BV31rmnOFPuElAUwcr3io09S60Lk0FeHVBU5/nrBM41jabfBLpggAABh1JREFUfj14NrW6KUuf7lff7lWf7k1N+vocriw/4P19ALep6qvKntezVwghhBBy5jPcOl9r9SvygyJyEsCTALxTRG5Y/fzbROT6VbUnA/jXAL5PRD68+np6UdvDX+ESQgjZPjaU+EJV3wbgbfv8/LMAnr76/n8BqKwQOplwReQKAL8JYAzgtap6nSs/C8AbATwewD0AnqWqn2xlMAPQLFtBn1V1n8dGCAGwnM36qpTr0rpSFpExgFcDeBqAiwE8R0QudtWeD+CLqvpwAL8O4BVtj4sQQkiP2ZBSbpMuPsO9BMAJVb1DVfcAvBnAla7OlQDesPr+LQCeKtFTx4QQQgbP0CbcLpTyeQA+bY5PAnhCXh1VnYnIlwF8C4C7bSURuRrA1avDU3Lb3R9rZcTbzyG4e0e+Ae9NPrw3+fDe5PPtbTR6Jz54w69ADtU8vZe/q60KmlplCjkCACJyLHqw+SDDe5MP700+vDf58N7k4zI6NYaqXtFGu5ukC6X8GQAXmOPzVz/bt46ITADcH8vgKUIIIWQQdDHh3gLgIhG5UER2ATwbwFFX5yiAq1bf/wiAv1TVnmWoIIQQQurTulJefSb7IgA3YPlY0OtU9VYReRmAY6p6FMuMHW8SkRMA/gHLSbmIVhJRDwTem3x4b/LhvcmH9yYf3puSCBeShBBCSPswtSMhhBDSAZxwCSGEkA7o/YQrIleIyMdF5ISIXLNP+Vki8ier8g+sNgM+EJS4N78gIsdF5KMicpOI/ONNjHMTFN0bU++HRURF5MA88lHm3ojIM1evnVtF5I+6HuOmKPE39VARebeIfGj1d1WYsH4IiMjrROQuEdk394Es+a3VffuoiDyu6zFuBara2y8sg6w+AeCfANgF8BEAF7s6LwTwmtX3zwbwJ5sed4/uzfcC+KbV9z/Ne5PUux+A9wK4GcDhTY+7L/cGwEUAPgTggavjb930uHt0b44A+OnV9xcD+OSmx93RvfkeAI8D8LGc8qcD+AssUyA/EcAHNj3mPn71fYXLtJD5FN4bVX23qt67OrwZy2egDwJlXjcA8HIs83Z/vcvBbZgy9+anALxaVb8IAKp6V8dj3BRl7o0C+ObV9/cH8NkOx7cxVPW9WD5BkseVAN6oS24G8AARObeb0W0PfZ9w90sLeV5eHVWdATidFnLolLk3ludj+Q70IFB4b1bK6wJVfWeXA+sBZV43jwDwCBF5n4jcvNrt6yBQ5t78CoDnrvZLvR7Az3QztN5T9f+jA8lWpXYk9RCR5wI4DOApmx5LHxCREYBXAXjehofSVyZYauVLsbQi7xWRR6vqlzY6qn7wHACvV9VfE5EnYZk/4FGqyj0fSSF9X+EyLWQ+Ze4NROQyANcCeIaqnupobJum6N7cD8CjALxHRD6J5WdORw9I4FSZ181JAEdVdaqqfwfgb7GcgIdOmXvzfAB/CgCq+n4AZ2O5scFBp9T/Rwedvk+4TAuZT+G9EZHHAvhdLCfbg/I5HFBwb1T1y6p6SFUfpqoPw/Lz7WeoaitJ2HtGmb+pt2O5uoWIHMJSMd/R5SA3RJl78ykATwUAEXkklhPuFzodZT85CuDHV9HKTwTwZVW9c9OD6hu9VsraXlrIrafkvXklgPsC+LNVHNmnVPUZGxt0R5S8NweSkvfmBgDfLyLHAcwB/JKqDt4albw3LwbweyLy81gGUD3vILzBF5E/xvJN2KHV59e/DGAHAFT1NVh+nv10ACcA3AvgJzYz0n7D1I6EEEJIB/RdKRNCCCGDgBMuIYQQ0gGccAkhhJAO4IRLCCGEdAAnXEIIIaQDOOESQgghHcAJlxBCCOkATriEdMxqP9XLV9//RxH5r5seEyGkfXqdaYqQgfLLAF4mIt8K4LEABp/9ixDCTFOEbAQR+Sss025eqqpf2fR4CCHtQ6VMSMeIyKMBnAtgj5MtIQcHTriEdIiInAvgDwFcCeCrB2hzd0IOPJxwCekIEfkmAG8F8GJVvQ3Ay7H8PJcQcgDgZ7iEEEJIB3CFSwghhHQAJ1xCCCGkAzjhEkIIIR3ACZcQQgjpAE64hBBCSAdwwiWEEEI6gBMuIYQQ0gH/H8Bw9yMUigYfAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Training Data for one function\n",
        "\n",
        "#### IC/BCC\n",
        "\n",
        "Ideally $G_\\theta(u)(\\textbf{y}_u)\\approx G(u)(\\textbf{y}_u)=0$, so compute the error → $\\mathcal{L}_{Operator}(\\theta)=\\frac{1}{NP}\\sum_{i=1}^N\\sum_{j=1}^P\\left|G_{\\theta}(u^{(i)})y_{uj}^{(i)}\\right|^2$"
      ],
      "metadata": {
        "id": "f6_C-hv7zG5p"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Sample points from the boundary and the inital conditions\n",
        "# Geneate subkeys \n",
        "subkeys = random.split(key, 4)\n",
        "# Here we regard the initial condition as a special type of boundary conditions\n",
        "x_bc1 = np.zeros((P_train // 3, 1))\n",
        "x_bc2 = np.ones((P_train // 3, 1))\n",
        "x_bc3 = random.uniform(key = subkeys[0], shape = (P_train // 3, 1))\n",
        "x_bcs = np.vstack((x_bc1, x_bc2, x_bc3))\n",
        "#Time\n",
        "t_bc1 = random.uniform(key = subkeys[1], shape = (P_train//3 * 2, 1)) # 200 points\n",
        "t_bc2 = np.zeros((P_train//3, 1)) # 100 points\n",
        "t_bcs = np.vstack([t_bc1, t_bc2])\n",
        "# Training data for BC and IC\n",
        "u_train = np.tile(u, (P_train,1)) # Add dimentions-> copy u P times\n",
        "y_train = np.hstack([x_bcs, t_bcs]) # stack the\n",
        "s_train = np.zeros((P_train, 1)) # Remember that the initial conditions are 0"
      ],
      "metadata": {
        "id": "ee4RwxTdx8Uw"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Collocation Points\n",
        "\n",
        "Select $Q$ points (i.e., collocation points) inside our domain → $\\textbf{y}_r=y_{r1},y_{r2},...,y_{rQ}$\n",
        "\n",
        "Compute the error related to the underlying governing laws (i.e., PDE)→$\\mathcal{L}_{Physics}(\\theta)=\\frac{1}{NQ}\\sum_{i=1}^{N}\\sum_{j=1}^{Q}\\left|R_{\\theta}^{(i)}(y_{r,j}^{(i)})-u^{(i)}(x_{r,j}^{(i)})\\right|^2$\n",
        "\n",
        "$$R_{\\theta}^{(i)}(y_{r,j}^{(i)})=\\frac{\\partial G_{\\theta}(u^{(i)})(y_{r,j}^{(i)})}{\\partial t}-D\\frac{\\partial^2 G_{\\theta}(u^{(i)})(y_{r,j}^{(i)})}{\\partial x^2}-k[G_{\\theta}(u^{(i)})(y_{r,j}^{(i)})]^2 $$"
      ],
      "metadata": {
        "id": "AyUVBQkP6KYN"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Sample collocation points\n",
        "x_r_idx= random.choice(subkeys[2], np.arange(Nx), shape = (Q_train,1),replace=False)\n",
        "x_r = x[x_r_idx]\n",
        "t_r = random.uniform(subkeys[3], minval = 0, maxval = 1, shape = (Q_train,1))\n",
        "\n",
        "# Training data for the PDE residual\n",
        "u_r_train = np.tile(u, (Q_train,1))\n",
        "y_r_train = np.hstack([x_r, t_r])\n",
        "f_r_train = u[x_r_idx] "
      ],
      "metadata": {
        "id": "VvPL-CyrzCal"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Training Data\n"
      ],
      "metadata": {
        "id": "BS3cApE9674l"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "\n",
        "We need to generate 5000 functions, so lets create a function to help us with the training and testing data\n"
      ],
      "metadata": {
        "id": "ORw11jP8xpgg"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Geneate training data corresponding to one input sample\n",
        "def generate_one_training_data(key, P, Q):\n",
        "    # Numerical solution\n",
        "    (x, t, UU), (u, y, s) = solve_ADR(key, Nx , Nt, P, length_scale)\n",
        "\n",
        "    # Geneate subkeys\n",
        "    subkeys = random.split(key, 4)\n",
        "\n",
        "    # Sample points from the boundary and the inital conditions\n",
        "    # Here we regard the initial condition as a special type of boundary conditions\n",
        "    x_bc1 = np.zeros((P // 3, 1))\n",
        "    x_bc2 = np.ones((P // 3, 1))\n",
        "    x_bc3 = random.uniform(key = subkeys[0], shape = (P // 3, 1))\n",
        "    x_bcs = np.vstack((x_bc1, x_bc2, x_bc3))\n",
        "\n",
        "    t_bc1 = random.uniform(key = subkeys[1], shape = (P//3 * 2, 1))\n",
        "    t_bc2 = np.zeros((P//3, 1))\n",
        "    t_bcs = np.vstack([t_bc1, t_bc2])\n",
        "\n",
        "    # Training data for BC and IC\n",
        "    u_train = np.tile(u, (P,1))\n",
        "    y_train = np.hstack([x_bcs, t_bcs])\n",
        "    s_train = np.zeros((P, 1))\n",
        "\n",
        "    # Sample collocation points\n",
        "    x_r_idx= random.choice(subkeys[2], np.arange(Nx), shape = (Q,1))\n",
        "    x_r = x[x_r_idx]\n",
        "    t_r = random.uniform(subkeys[3], minval = 0, maxval = 1, shape = (Q,1))\n",
        "\n",
        "    # Training data for the PDE residual\n",
        "    '''For the operator'''\n",
        "    u_r_train = np.tile(u, (Q,1))\n",
        "    y_r_train = np.hstack([x_r, t_r])\n",
        "    '''For the function'''\n",
        "    f_r_train = u[x_r_idx]\n",
        "    return u_train, y_train, s_train, u_r_train, y_r_train, f_r_train"
      ],
      "metadata": {
        "id": "ffEaJ96I704M"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# GRF length scale\n",
        "length_scale = 0.2\n",
        "\n",
        "# Resolution of the solution\n",
        "Nx = 100\n",
        "Nt = 100\n",
        "\n",
        "N = 5000 # number of input samples\n",
        "m = Nx   # number of input sensors\n",
        "P_train = 300 # number of output sensors, 100 for each side \n",
        "Q_train = 100  # number of collocation points for each input sample"
      ],
      "metadata": {
        "id": "hHB7XD_NASf9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Select N keys to create N Functions\n",
        "key = random.PRNGKey(0)\n",
        "keys = random.split(key, N)"
      ],
      "metadata": {
        "id": "tPz8y43RAiBt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "config.update(\"jax_enable_x64\", True)\n",
        "u_train, y_train, s_train, u_r_train, y_r_train, f_r_train = vmap(generate_one_training_data, (0, None, None))(keys, P_train, Q_train)"
      ],
      "metadata": {
        "id": "hROe7KViAZ5E"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Reshape the data \n",
        "u_bcs_train = np.float32(u_train.reshape(N * P_train,-1))\n",
        "y_bcs_train = np.float32(y_train.reshape(N * P_train,-1))\n",
        "s_bcs_train = np.float32(s_train.reshape(N * P_train,-1))\n",
        "\n",
        "u_res_train = np.float32(u_r_train.reshape(N * Q_train,-1))\n",
        "y_res_train = np.float32(y_r_train.reshape(N * Q_train,-1))\n",
        "f_res_train = np.float32(f_r_train.reshape(N * Q_train,-1))"
      ],
      "metadata": {
        "id": "EiWKTB66AxPg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Testing Data"
      ],
      "metadata": {
        "id": "6Yv2UV8dGTg0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Geneate test data corresponding to one input sample\n",
        "def generate_one_test_data(key, P):\n",
        "    Nx = P\n",
        "    Nt = P\n",
        "    (x, t, UU), (u, y, s) = solve_ADR(key, Nx , Nt, P, length_scale)\n",
        "\n",
        "    XX, TT = np.meshgrid(x, t)\n",
        "\n",
        "    u_test = np.tile(u, (P**2,1))\n",
        "    y_test = np.hstack([XX.flatten()[:,None], TT.flatten()[:,None]])\n",
        "    s_test = UU.T.flatten()\n",
        "\n",
        "    return u_test, y_test, s_test"
      ],
      "metadata": {
        "id": "lZ0Pc5m2_xP2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "N_test = 100 # number of input samples \n",
        "key = random.PRNGKey(12345) # it should be different than the key we used for the training data\n",
        "P_test = 100\n",
        "Nx = m\n",
        "keys = random.split(key, N_test)\n",
        "\n",
        "config.update(\"jax_enable_x64\", True)\n",
        "u_test, y_test, s_test = vmap(generate_one_test_data, (0, None))(keys, P_test)\n",
        "\n",
        "#Reshape Data\n",
        "u_test = np.float32(u_test.reshape(N_test * P_test**2,-1))\n",
        "y_test = np.float32(y_test.reshape(N_test * P_test**2,-1))\n",
        "s_test = np.float32(s_test.reshape(N_test * P_test**2,-1))"
      ],
      "metadata": {
        "id": "v6Ng55PJGq2F"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "config.update(\"jax_enable_x64\", False)"
      ],
      "metadata": {
        "id": "QGV6BUeKh_ai"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Data generator\n",
        "class DataGenerator(data.Dataset):\n",
        "    def __init__(self, u, y, s, \n",
        "                 batch_size=64, rng_key=random.PRNGKey(1234)):\n",
        "        'Initialization'\n",
        "        self.u = u # input sample\n",
        "        self.y = y # location\n",
        "        self.s = s # labeled data evulated at y (solution measurements, BC/IC conditions, etc.)\n",
        "        \n",
        "        self.N = u.shape[0]\n",
        "        self.batch_size = batch_size\n",
        "        self.key = rng_key\n",
        "\n",
        "    def __getitem__(self, index):\n",
        "        'Generate one batch of data'\n",
        "        self.key, subkey = random.split(self.key)\n",
        "        inputs, outputs = self.__data_generation(subkey)\n",
        "        return inputs, outputs\n",
        "\n",
        "    @partial(jit, static_argnums=(0,))\n",
        "    def __data_generation(self, key):\n",
        "        'Generates data containing batch_size samples'\n",
        "        idx = random.choice(key, self.N, (self.batch_size,), replace=False)\n",
        "        s = self.s[idx,:]\n",
        "        y = self.y[idx,:]\n",
        "        u = self.u[idx,:]\n",
        "        # Construct batch\n",
        "        inputs = (u, y)\n",
        "        outputs = s\n",
        "        return inputs, outputs"
      ],
      "metadata": {
        "id": "Jz0NWj_QMcZe"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#DeepONet"
      ],
      "metadata": {
        "id": "5kcTUiVJMfhL"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Define the neural net\n",
        "def MLP(layers, activation=relu):\n",
        "  ''' Vanilla MLP'''\n",
        "  def init(rng_key):\n",
        "      def init_layer(key, d_in, d_out):\n",
        "          k1, k2 = random.split(key)\n",
        "          glorot_stddev = 1. / np.sqrt((d_in + d_out) / 2.)\n",
        "          W = glorot_stddev * random.normal(k1, (d_in, d_out))\n",
        "          b = np.zeros(d_out)\n",
        "          return W, b\n",
        "      key, *keys = random.split(rng_key, len(layers))\n",
        "      params = list(map(init_layer, keys, layers[:-1], layers[1:]))\n",
        "      return params\n",
        "  def apply(params, inputs):\n",
        "      for W, b in params[:-1]:\n",
        "          outputs = np.dot(inputs, W) + b\n",
        "          inputs = activation(outputs)\n",
        "      W, b = params[-1]\n",
        "      outputs = np.dot(inputs, W) + b\n",
        "      return outputs\n",
        "  return init, apply"
      ],
      "metadata": {
        "id": "SiTZDAclMjfX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**In summary:**\n",
        "\n",
        "*   Let $y=(x,t)$.\n",
        "*   Also remember that $u$ is only a function of $x$.\n",
        "\n",
        "To train a PI-DeepOnet:\n",
        "\n",
        "1.   Select $N$ source terms (i.e, input functions) →$\\{u^{(i)}(\\textbf{x})\\}_{i=1}^{N}={[u^{(1)}(\\textbf{x}),u^{(2)}(\\textbf{x}),...,u^{(N)}(\\textbf{x})]}$.\n",
        "2.   Evaluate our $N$ functions at $m$ points (i.e., input sensors) →$\\begin{bmatrix}u^{(1)}(x_1)&u^{(1)}(x_2)&...&u^{(1)}(x_m)\\\\\n",
        "u^{(2)}(x_1)&u^{(2)}(x_2)&...&u^{(2)}(x_m)\\\\\n",
        "\\vdots&\\vdots&\\ddots&\\vdots\\\\\n",
        "u^{(N)}(x_1)&u^{(N)}(x_2)&...&u^{(N)}(x_m)\\end{bmatrix}$\n",
        "3.   Send the $m$ outputs of our $N$ functions to our **branch network** → $b_k\\begin{pmatrix}u^{(1)}(x_1)&u^{(1)}(x_2)&...&u^{(1)}(x_m)\\\\\n",
        "u^{(2)}(x_1)&u^{(2)}(x_2)&...&u^{(2)}(x_m)\\\\\n",
        "\\vdots&\\vdots&\\ddots&\\vdots\\\\\n",
        "u^{(N)}(x_1)&u^{(N)}(x_2)&...&u^{(N)}(x_m)\\end{pmatrix}$\n",
        "4.   Select $P$ points (i.e., output sensors) from our boundary and initial conditions → $\\textbf{y}_u=y_{u1},y_{u2},...,y_{uP}$\n",
        "5.  Send our output sensors to our **trunk network**→$t_k(y_{u1},y_{u2},...,y_{uP})$\n",
        "6. Approximate our operator by computing the dot product between the output of our **branch network** and the output of our **trunk network**→ $G_\\theta(u)(\\textbf{y})=\\sum_{k=1}^q\\underset{Branch}{\\underbrace{b_k\\left(u(x_1),u(x_2),...,u(x_m)\\right)}}.\\underset{Trunk}{\\underbrace{t_k(\\textbf{y})}}$\n",
        "7. Ideally $G_\\theta(u)(\\textbf{y}_u)\\approx G(u)(\\textbf{y}_u)=0$, so compute the error → $\\mathcal{L}_{Operator}(\\theta)=\\frac{1}{NP}\\sum_{i=1}^N\\sum_{j=1}^P\\left|G_{\\theta}(u^{(i)})y_{uj}^{(i)}\\right|^2$\n",
        "8. Select $Q$ points (i.e., collocation points) inside our domain → $\\textbf{y}_r=y_{r1},y_{r2},...,y_{rQ}$\n",
        "9. Compute the error related to the underlying governing laws (i.e., PDE)→$\\mathcal{L}_{Physics}(\\theta)=\\frac{1}{NQ}\\sum_{i=1}^{N}\\sum_{j=1}^{Q}\\left|R_{\\theta}^{(i)}(y_{r,j}^{(i)})-u^{(i)}(x_{r,j}^{(i)})\\right|^2$\n",
        "\n",
        "10. Compute the total loss→$\\mathcal{L}(\\theta)=\\mathcal{L}_{operator}(\\theta)+\\mathcal{L}_{Physics}(\\theta)$\n",
        "11. Update our NN parameters (i.e., branch and trunk) to minimize  $\\mathcal{L}(\\theta)$.\n",
        "12. Repeat the process until $G_\\theta(u)(x,t)\\approx G(u)(x,t)$.\n",
        "\n",
        "The residual will be:\n",
        "\n",
        "$$R_{\\theta}^{(i)}(y_{r,j}^{(i)})=\\frac{\\partial G_{\\theta}(u^{(i)})(y_{r,j}^{(i)})}{\\partial t}-D\\frac{\\partial^2 G_{\\theta}(u^{(i)})(y_{r,j}^{(i)})}{\\partial x^2}-k[G_{\\theta}(u^{(i)})(y_{r,j}^{(i)})]^2 $$"
      ],
      "metadata": {
        "id": "oxp4n-iJ8YG1"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Define the model\n",
        "class PI_DeepONet:\n",
        "    def __init__(self, branch_layers, trunk_layers):    \n",
        "        # Network initialization and evaluation functions\n",
        "        self.branch_init, self.branch_apply = MLP(branch_layers, activation=np.tanh)\n",
        "        self.trunk_init, self.trunk_apply = MLP(trunk_layers, activation=np.tanh)\n",
        "\n",
        "        # Initialize\n",
        "        branch_params = self.branch_init(rng_key = random.PRNGKey(1234))\n",
        "        trunk_params = self.trunk_init(rng_key = random.PRNGKey(4321))\n",
        "        params = (branch_params, trunk_params)\n",
        "\n",
        "        # Use optimizers to set optimizer initialization and update functions\n",
        "        self.opt_init, \\\n",
        "        self.opt_update, \\\n",
        "        self.get_params = optimizers.adam(optimizers.exponential_decay(1e-3, \n",
        "                                                                      decay_steps=2000, \n",
        "                                                                      decay_rate=0.9))\n",
        "        self.opt_state = self.opt_init(params)\n",
        "\n",
        "        # Used to restore the trained model parameters\n",
        "        _, self.unravel_params = ravel_pytree(params)\n",
        "\n",
        "        self.itercount = itertools.count()\n",
        "\n",
        "        # Loggers\n",
        "        self.loss_log = []\n",
        "        self.loss_bcs_log = []\n",
        "        self.loss_res_log = []\n",
        "\n",
        "    # Define DeepONet architecture\n",
        "    def operator_net(self, params, u, x, t):\n",
        "        branch_params, trunk_params = params\n",
        "        y = np.stack([x, t])\n",
        "        B = self.branch_apply(branch_params, u)\n",
        "        T = self.trunk_apply(trunk_params, y)\n",
        "        outputs = np.sum(B * T)\n",
        "        return  outputs\n",
        "  \n",
        "    # Define ODE/PDE residual\n",
        "    def residual_net(self, params, u, x, t):\n",
        "        s = self.operator_net(params, u, x, t)\n",
        "        s_t = grad(self.operator_net, argnums=3)(params, u, x, t)\n",
        "        s_x = grad(self.operator_net, argnums=2)(params, u, x, t)\n",
        "        s_xx= grad(grad(self.operator_net, argnums=2), argnums=2)(params, u, x, t)\n",
        "\n",
        "        res = s_t - 0.01 * s_xx - 0.01 * s**2 \n",
        "        return res\n",
        "\n",
        "    # Define boundary loss\n",
        "    def loss_bcs(self, params, batch):\n",
        "        inputs, outputs = batch\n",
        "        u, y = inputs\n",
        "\n",
        "        # Compute forward pass\n",
        "        s_pred = vmap(self.operator_net, (None, 0, 0, 0))(params, u, y[:,0], y[:,1])\n",
        "\n",
        "        # Compute loss\n",
        "        loss = np.mean((outputs.flatten() - s_pred)**2)\n",
        "        return loss\n",
        "\n",
        "    # Define residual loss\n",
        "    def loss_res(self, params, batch):\n",
        "        inputs, outputs = batch\n",
        "        u, y = inputs\n",
        "        # Compute forward pass\n",
        "        pred = vmap(self.residual_net, (None, 0, 0, 0))(params, u, y[:,0], y[:,1])\n",
        "\n",
        "        # Compute loss\n",
        "        loss = np.mean((outputs.flatten() - pred)**2)\n",
        "        return loss   \n",
        "\n",
        "    # Define total loss\n",
        "    def loss(self, params, bcs_batch, res_batch):\n",
        "        loss_bcs = self.loss_bcs(params, bcs_batch)\n",
        "        loss_res = self.loss_res(params, res_batch)\n",
        "        loss = loss_bcs + loss_res\n",
        "        return loss \n",
        "\n",
        "    # Define a compiled update step\n",
        "    @partial(jit, static_argnums=(0,))\n",
        "    def step(self, i, opt_state, bcs_batch, res_batch):\n",
        "        params = self.get_params(opt_state)\n",
        "        g = grad(self.loss)(params, bcs_batch, res_batch)\n",
        "        return self.opt_update(i, g, opt_state)\n",
        "\n",
        "    # Optimize parameters in a loop\n",
        "    def train(self, bcs_dataset, res_dataset, nIter = 10000):\n",
        "        # Define data iterators\n",
        "        bcs_data = iter(bcs_dataset)\n",
        "        res_data = iter(res_dataset)\n",
        "\n",
        "        pbar = trange(nIter)\n",
        "        # Main training loop\n",
        "        for it in pbar:\n",
        "            # Fetch data\n",
        "            bcs_batch= next(bcs_data)\n",
        "            res_batch = next(res_data)\n",
        "\n",
        "            self.opt_state = self.step(next(self.itercount), self.opt_state, bcs_batch, res_batch)\n",
        "            \n",
        "            if it % 100 == 0:\n",
        "                params = self.get_params(self.opt_state)\n",
        "\n",
        "                # Compute losses\n",
        "                loss_value = self.loss(params, bcs_batch, res_batch)\n",
        "                loss_bcs_value = self.loss_bcs(params, bcs_batch)\n",
        "                loss_res_value = self.loss_res(params, res_batch)\n",
        "\n",
        "                # Store losses\n",
        "                self.loss_log.append(loss_value)\n",
        "                self.loss_bcs_log.append(loss_bcs_value)\n",
        "                self.loss_res_log.append(loss_res_value)\n",
        "\n",
        "                # Print losses\n",
        "                pbar.set_postfix({'Loss': loss_value, \n",
        "                                  'loss_bcs' : loss_bcs_value, \n",
        "                                  'loss_physics': loss_res_value})\n",
        "           \n",
        "    # Evaluates predictions at test points  \n",
        "    @partial(jit, static_argnums=(0,))\n",
        "    def predict_s(self, params, U_star, Y_star):\n",
        "        s_pred = vmap(self.operator_net, (None, 0, 0, 0))(params, U_star, Y_star[:,0], Y_star[:,1])\n",
        "        return s_pred\n",
        "\n",
        "    @partial(jit, static_argnums=(0,))\n",
        "    def predict_res(self, params, U_star, Y_star):\n",
        "        r_pred = vmap(self.residual_net, (None, 0, 0, 0))(params, U_star, Y_star[:,0], Y_star[:,1])\n",
        "        return r_pred"
      ],
      "metadata": {
        "id": "AYvrwHItMlXJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Initialize DeepONet"
      ],
      "metadata": {
        "id": "gPuIzH2_NB5T"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Initialize model\n",
        "branch_layers = [m, 50, 50, 50, 50, 50]\n",
        "trunk_layers =  [2, 50, 50, 50, 50, 50]\n",
        "model = PI_DeepONet(branch_layers, trunk_layers)"
      ],
      "metadata": {
        "id": "jblmIk1kNmzX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Create data set\n",
        "batch_size = 10000\n",
        "bcs_dataset = DataGenerator(u_bcs_train, y_bcs_train, s_bcs_train, batch_size)\n",
        "res_dataset = DataGenerator(u_res_train, y_res_train, f_res_train, batch_size)"
      ],
      "metadata": {
        "id": "heAqQrjQNSTE"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Train"
      ],
      "metadata": {
        "id": "xPy_LLl4N7NW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Train\n",
        "model.train(bcs_dataset, res_dataset, nIter=40000)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "edzeJzhZN1zB",
        "outputId": "7b52adc8-b0ad-4d4c-e5bb-2a7dfbe50c4a"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "100%|██████████| 120000/120000 [45:38<00:00, 43.82it/s, Loss=6.180332e-05, loss_bcs=1.8806511e-05, loss_physics=4.2996813e-05]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Evaluate"
      ],
      "metadata": {
        "id": "mHE436wq13E-"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Results for 100 test functions"
      ],
      "metadata": {
        "id": "9se1Wt2_8OKa"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "params = model.get_params(model.opt_state)\n",
        "s_pred = model.predict_s(params, u_test, y_test)[:,None]\n",
        "error_s = np.linalg.norm(s_test - s_pred) / np.linalg.norm(s_test)\n",
        "print(error_s)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "db01GMpB32SR",
        "outputId": "297f48c6-360c-4d9f-f247-f51fec107dd5"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0.0071554612\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Result for 1 test function"
      ],
      "metadata": {
        "id": "Jy0GuJ7h8RUe"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "N_test = 1 # number of input samples \n",
        "key = random.PRNGKey(456)\n",
        "P_test = 100\n",
        "Nx = m\n",
        "keys = random.split(key, N_test)\n",
        "\n",
        "config.update(\"jax_enable_x64\", True)\n",
        "u_test, y_test, s_test = vmap(generate_one_test_data, (0, None))(keys, P_test)\n",
        "\n",
        "#Reshape Data\n",
        "u_test = np.float32(u_test.reshape(N_test * P_test**2,-1))\n",
        "y_test = np.float32(y_test.reshape(N_test * P_test**2,-1))\n",
        "s_test = np.float32(s_test.reshape(N_test * P_test**2,-1))"
      ],
      "metadata": {
        "id": "w-9at0SY2Z27"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Predict\n",
        "params = model.get_params(model.opt_state)\n",
        "s_pred = model.predict_s(params, u_test, y_test)\n"
      ],
      "metadata": {
        "id": "qUftmA4i14QC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Generate an uniform mesh\n",
        "x = np.linspace(0, 1, Nx)\n",
        "t = np.linspace(0, 1, Nt)\n",
        "XX, TT = np.meshgrid(x, t)"
      ],
      "metadata": {
        "id": "TlJZflzE1_cM"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Grid data\n",
        "S_pred = griddata(y_test, s_pred.flatten(), (XX,TT), method='cubic')\n",
        "S_test = griddata(y_test, s_test.flatten(), (XX,TT), method='cubic')"
      ],
      "metadata": {
        "id": "p8F5QSUl2Bhz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Real\n",
        "plot(XX,TT, S_test)\n",
        "#Prediction\n",
        "plot(XX,TT, S_pred)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 721
        },
        "id": "_RGIIFqW3XDX",
        "outputId": "2d8ad339-c9b5-41bd-af66-68389b619efd"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 504x360 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFgCAYAAADgoJN2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO29fax151need1773PsfDgJ4SUEOYakwmnjSaABTwJiREJJGCej2n8wBSfKTNMGPC2YYQggBYECCvPBRwullQd4oVEADRiIKvQijNyKhkYTkYxfFJpi0zCuoYkDU3/guKEm77s/7vljnzfnfq7nrHuttb/O2udcP+lIe51nfe2119rPXte6nus2d4cQQgghVmd02jsghBBC7DvqTIUQQog1UWcqhBBCrIk6UyGEEGJN1JkKIYQQa6LOVAghhFgTdaZCCCHEmqgzFUIIIdZEnakQG8DMvsDM/pWZPWVm7zWz/8PM/pcNb+P/MbP/apPrFEJsBlMCkhDrY2Y/AeB6d/82M/sCAH8A4Evd/a/WXO+jAP62u3/UzL4JwDe7+zduYJeFEBtEd6ZCbIY3APj1o9dvB3DfBjrSCwC+EMBDR/+6BODrzOzF66xXCLF51JkKsQZmdmhmTwN4FYDfNLN/B+BNAP5NmOfHzOw3wvSPm9nvmNlhst4vBfBJLK/RJ83sSQAzAL8P4L/dzrsRQqzK5LR3QIh9xt2vmtlXA/iAu38hAJjZ4wA+Hmb7UQCPmNmrAbwWwG0A/ht3v5qs92Ez+x4AX+vu33zt/2b2RwC+fAtvRQixBupMhVifvwng34bpFwD4zLUJd3/SzH4SwC8AeD6WHenTHdb75Vg+e418BsAXrbe7QohNI5lXiPXhzvQpADfQPB/FUgr+Pnf/5IrrxdF6P73KTgohtoc6UyHW58tRdnofA/DyaxNm9ioAP43lnenf77JCMxsBeCXqO9NXoO5ghRCnjDpTIdaHO9P7ALwOAMzsRgC/CeAfAPg2AK8ys9fHhc3sfWb2Plrns47+RmG+6wF8JYB/tdndF0KsizpTIdbgaJjK5wH49+HfvwjgzWb2fCw71p9w90vu/gyAHwfwv9FqbgLwofgPd/8vAH4GwENHY00B4G8D+F13/7PNvxMhxDootEGILWBm/zuAx9z9n7TMd4jlXe2Xufu0Zd6PAHiHu//h5vZUCLEJ1JkKIYQQa7J1mfcop/QxMzvx17Qt+adm9rCZfczMvmLb+ySEEEJskl08M30floPUm3gTgJuP/u7C0vUohBBCbAUzu83MPn50E/euE9q/2Mw+YGYfPbrJe3PbOrfembr7BwH8RTLLHQB+0Zd8GMALzEyD0oUQQmwcMxsDuAfLG7lbALzFzG6h2X4AwK+5+6sB3Ang/2xb7xASkG7EMoP0Go8e/e/PeUYzuwvLu1cc4Dlf+cLx3zhqKOdbhJ8IizG1hXc8H5fPi+cHoOnj9jkdqdkkLDsp1zOm9U7Gi+O2Udk2CtNjozaaNnjSFpZD83NwS9rWwflDWJFs/7It8HLFsaK2kSfzevO848WibFuENqe2eTk9DtNGbYjT3DZL5q3avLkt7Gu1jQUd89i8jqcifmBGn97ITn4NAJNwAY+TNgA4GDe3TcYnzwdgTvPOwryzcTlvnJ5ZudzUQhvdm0yd1hO+mGZz2v78+H3OZuV79lk57yS0j6flvOPZya+X08fzjuZlG09bOAeMTwGa/k/z33/C3b8AG+Y2M39ixWV/H7jf3ZsU0dcAeNjdHwEAM7sXy5u6h8I8DuB5R6+fD6DVQT+EzrQz7n4RwEUAePHkVv8fXnAZQN1hXn3W8etnXlB+8n/5Qg+vyy+Vp15cTn/6Rcdn2Ke/oDzbnnpRMF6+sDRhvuD55fTnPe84gvV5zynjWG+4/nj6hsOy7VmT8mq4Plwd14/KtuvsePrQyn2NnckE9EW6IeYb6kzHWSdI0+PwXiYo3/N1fjx93aL8PJ5F07H9OdPyM3h2mH7Olc8WbTf81fH0c/+qbHveZ54ppp/7meMCModPl214Kkw/TYVmnvwvNG9o/wua9+nPNrf95ZXj15+hSOC/pOkr4VhepW9Z7ogzYgd2HV2kzwq/XJ9Lef8vuP749QufVba98Nnl9Bc+9/j1hedQWwihetHziqb/TNNPvPB4+vEXlG2P33A8/cR1zy3a/tPB8TYe87Lt8Wk5/dgzx/v++F+W7+vJT1933PbkdUXb4ony+HzeY8fH7vP/v/Ir/IYnjzve5z9edsLP/Yvj6Wc/XV6vz6Zwy8NnjtsndHpwx/uPnrT/iC3wBIDLKy5rwIWk+aQbuNfSPD8E4F+a2XcAeA6WVaFShjDO9FNYjrO7xkuO/ieEEOI8Mx6t9gdcMLPL4e+unlt+C4D3uftLALwZwC8dpZI1MoQ700sA7j661X4tgKfdvZJ4hRBCnCMMtczflTmecPdbG1q73MC9A0fGWXf/vaP0sQsAHmva5NY7UzP7FQCvx/KXwqMAfhDAwdFO/gyWCTFvBvAwgGcA/L1t7g9LFKsSldS2p0rl4yqjNjvx9XK9zSeS07yL8KNpTs/v4mr4uc6mnqH2eWa6CPOydDsLbbxvByRRj8K8C9p+nJ7TD8qZlZLjJDwgmo3KtkV41jfPnqWN6XkZzTsNkufhhCTPg+Zne+n0IbXF6WfRpR3l2kO6CJ5FZoGMeY8vt+s67itLwLHtcNLcBly7C1mSHauD8vNZ0HPaOD0b0bzhXOJzZ4bjab4G5t48zZfoInwvLOg7gp7aYJx8h2VtkU19D24XKz/fXqRv8AEAN5vZy7DsRO8E8Faa5xMAvh7A+8zsFQCuB/B4ttKtd6bu/paWdgfw7dveDyGEEHuEAZhsxosRcfeZmd0N4H4AYwDvdfcHzew9AC67+yUA3w3g58zsu7C8X3q7tyQcDUHmFUIIIXaGu9+HpSoa//fu8PohAF/TZ51nvjPtI2ewRDIOklYmnzhJXyzhZES5tpIqWfaN81oiF/MwhPh7ikcodN/Vcns9ZF2Wu4o2Wk8cHsT7Nu2xt1OLDmaS7fjYWfOxi3JtNZwhSFALkgZnNPRingzTSGXeTB5dtY1l1Tn94M4cu31k3q5ybdq2xvEIx3JB0vqMpuPn7CzzJudHlHZnNBSGZd/skU58/MPfJxOSfUfh48m+39qGvwwewxoy7+45852pEEKIPWVVA9IpoM5UCCHE8LB1DEi7R53pioyj+47aKhffolneSZdjx250qPK8QX6aVM7WAD9CX9V5nryP5TZXXHEPSXrR0c1bJdOQu3cSHJpTkvgOw7ws8c3DvOzenZOMWLhFqW3SVQIGesi8LJ3OmtsymZe/zK5SrE6E543TvK9dnb5tMm8mkScyPH9esZ3dvNHhzY8limuyz2MadvWHaVbZr+vx2GjUR4YfOpJ5hRBCiA0gmVcIIYRYA8m8+0ObJNLV/dYWVxolnGqwdhzI3SadJo7d0ulbnoCjIPRWu7piZkObjLtq8H2UyhYcOs87a81t8zCYfkHvek4D7xdhgHfVFo4zBzpEOXA+Zqdvs3t0SnLkZNXQhj4ScAxxyIL2mStscQ8HneVhvouIxyTL5uWAiShDZ7Iut3Nbcsw5tCF+fpWcH4M76NqK1wFfv9VjmjDNbX1GAMTvrVGyHH+/xe+zyulL6n02707RnakQQgixBoa6EtCA2Z89FUIIIQbKubszzQc5N0sKLKfE6cyZ1waXkszaojTEaaqFs5Xl4iBNjUjy7FOQLZNuF8nvsszBzIyK0IZ8X+N75kCHKPtW7l3S1GL7IUl8UdqtHLtFNm/3UAB2+ha5sX3cvJyp+9npyfMBJJ3SRcDzZsx6yG7xroL3dRPhE0DnXGM+5jxd5CyznB+kf36EEs87fvRR53Afv/ZFsyTM3x/8vWSptNvctneYSeYVQggh1kYGJBHxjjFi/IuUDQ2emR2KKilZnGCzOamNVe8+exUOD/u6QHmHMKHyGVMPd9zG1WeO28a03zzu9CCajMiAFMedXl/dtSbjTNnEEu6EOGowRt2NsooyPF3dtca7T3KUpOakHlVj+twpZEaibcQJVscqjjOlz47HmUaTEc3r8fygc8eTsd/ptU6qU1k1pmyzROmyTE3b0l3qzu5+Nc5UCCGEWBPJvEIIIcQG0J3p6dJVhmirqlBWZ2g2BbAsk01nUm6bTFSMMyXp1DIJOE5WcYLbkW7bxsw2URQHJ+kW3iz7RsmXl52TNFdJsMWY1PIkiNL2jMcYZlGDiSGprjBz3HbIUuX1JMFeCfLtZ0nKjcumxcF7xAky8+TLLRtnuqrJqM+4Wy74XUjr+TjTwiBWGZCaTUZlxGfzGFSArvVMuqW2PtWsulaR2QujkkF3pkIIIcR67FcC0v7sqRBCCDFQzvyd6WnIGTyGLJJVW6lkXR63lriCR3GcKUmOo6DisfmP92eMKI8272tr9GGybHYMiuLgrD7SYqPwW9CqOMEozZWy5YxcwrPQzmNSD4tqIjx2NMYJNo9bXE43O39j1N3hAf2+7ePmjdLulMeSxqoxXBw8kXX5ziDO2ydOMB33mryPNcaZLhK3dVU1JhYAZyk3nBP8iCA+XsiuV4Ae01Rj049fjzcUc9rnu2+Q1WYk8wohhBBroqExQgghxAbQnelw6SN9xAHRHOEVHXfZoGogjxfMqsbUMlFz28Sa24rqKi1VYmYdAxbaJK1s3pT4Hmm5A5JriwhBloA9iROsYgpj4AU5bZPC4dH1WRWUpl/UUVbk0IYY6LAg12ka4tAnWq+IIWQ3b49QyT5y4Mpu3h5VY5L3HOXz6phn0ZCJm5cfEcTQhqwYOFA6eKtAluQ7YpS1VW7ezYQ4DMLtqxJsQgghxAZQZyqEEEKsgQxI+0NbaENX2H3Hsm4h7yRtTCUFRUk4GTxurHkW2byNm2tlVSm3T4DDPAlt4GCGKPtWsnOQ5rxy85IEG6vGoHngPWcTF/Jf5fRtlhGriiVJuEAV4rBycfAQ/sCybp9s3j4y76qO3VWLpSeVYZxk+CqrN7RXARzIrrvg9OXHMlXVmCzo5bitqlBVhcl0+wwGIdWug2ReIYQQYgPs0Z3p/nT7QgghxEDRnWlHWMrlTMxIZo6sSi8lQQwZWW7vqBp0HiZa3LzpNntIub0cvA2w85jVLZZ9i2WDRMzzTYzdvNGtmUjALBWGtuk4d4tGWZEDHaIkPD1gmTcJcegjAWcu4E1l8zKbKAC+hps3OnhZPq9LsIXPsvqck0CH5NELy76eZHTHQ94m41qRsUsBLVsIdDg1tjjO1MxuA/BTAMYAft7df4TafxLA1x1NPhvAi9z9Bdk61ZkKIYQYHlsqwWZmYwD3AHgjgEcBPGBml9z9oWvzuPt3hfm/A8Cr29YrmVcIIcQwGY9W+8t5DYCH3f0Rd78K4F4AdyTzvwXAr7StVHemCZmUG2H3HRMlncy9W0nA7AYMUiFLSHF6lGTzdjQCnkiUcttk3LR8W8dsXt7GiDVqa3bzRml3TMeDpdxJmK6yeaNj10ieDXLgdVV5NtpmEdrQLEey0zcNcegj5U7D9LSlBFtka9m8YZrLzPVx8yZO6MIlnYRoLKdjxi6X7AvrYUd3OM+mCzoHkkcx/D0Qr3Wj750sFIbDZDJWDXQ4NUl4e0NjbgTwyTD9KIDXnrgLZl8C4GUA/nXbStWZCiGEGCBrDY25YGaXw/RFd7+4wnruBPB+d2/9SaHOVAghxPBY7870CXe/taHtUwBuCtMvOfrfSdwJ4Nu7bPBMdKZ9JIpS6vCkrWU9mVxLbZ5IMbGNpZ9qEHiRA0rrSTJC+2TzZiyybRDzxOWYMQvLcVk1JrqWOce3aGN51sv1ToJEe5CU32IpOcq+VabvqFlWrAIdOkqVAIU4XEeX72ejdEptV+fNbaeSzduxzFqP0AZ2Qs8Lt3VzybXldJDa2c2bnGfZNTln2Tdc39V3RDiuE2rjx01ZibZR4vTdO7bn5n0AwM1m9jIsO9E7Aby12rzZ3wDweQB+r8tKz0RnKoQQ4qyxnQQkd5+Z2d0A7sdyaMx73f1BM3sPgMvufulo1jsB3OvunW5B1JkKIYQ4V7j7fQDuo/+9m6Z/qM86z3Vn2uZS6zwAOpFsgFLyzPI6q+CDxLE7yZyCmazaQ/np48jNpNw+YRQRik89wXTa/Ku1cAXzsUqyV+dUYmuG45OgcvoGObAq21VNBxmRfm3PEgk4DXFgyTO6Yqd08sYsXA6GnW8pmzeTeTcRPgEAB81O6OiSXnD5PJbaMzk/nGdZMEN9/aJ5Xv4IwjR/Gtloga4jDpi9CW3gL4EBc647UyGEEANGQfdCCCHEGmxvnOlWOJOd6TYkDJZyV91GJQVFh19Sco2nMwm4lj8XDfOVciiThStUpajS/aG9SR7nR1WnfV+P3xcHVcTQBg57GJG7d56UciuyeTmYIczLea5ZCTbOhS1LsJEcmYQ4jFaVR7kty+blO4OrMzTC844TSbprNm+LzLvoeOwqaX3E02E9fH6Ec4tznudhuroGKld/c+hKVoKN6Tp6YZ0Sk8OQgbdjQNoWZ7IzFUIIsefozlQIIYTYACPdme4lecBDiSVZmjOSadjBG8nl0EQSHrHMWuxd0cYSaLn9xqZWKbdcT+5kLNabruf4tVWyLi8XLjQ65qOQTlEN0Gc5Lqx3RlJuDHFYcDCDZdJgs6xY58IGObIltCFOH05WlHnbSrBFrvJFEY4df8jsvMzcvDE4Yg03b3TpsmN3lkjrdWiDNbbNgsO7KsEWzw9qa3P3NlGXVaPpRXPbqgxD1iW2VDVmW+xPty+EEEIMFN2ZCiGEGB7bixPcCueuM82zK5t1GJZyI2OScafUHsfI83j5SFWWKRn0zbJqDHFY0PvI5FlW5rpKuW2SVeY87gy7mVn2DdcZS9mxNJZRGhi7e2M2L5fYiiEOMcABKJ2+1/P2kxAHdpJGpynLkVmIw2FSjiyVR7nkWSbzMquGNvRx88bMYd7XSfPx4GMVJfGqjT6f6MbmcA5PQj3iuT1rKcE2T7J54/cCl1xbtcxa6vpNTNmDYo9k3nPXmQohhNgDzGRAEkIIIdZGd6ZnnywTM3Pv1uXZkiCGFXN8+bdcpuJlbmImSrttOcJFW5Lb2ydEolLhC/mrfNcW3LxVjjE7OcOyWS5rlc0b5WEOdOAQiejYZckxuk6TfFkAmIbSZbMrpVY3yaTcmNXLGj2XZMvo8+U2Sty8WSm52JY5lsHl65qDM7KSeEBZpo9d7PEcyMJJ+Nqurt9C5i2aCmm3reSadZRymUE6djP27JnpTvbUzG4zs4+b2cNm9q4T2r/YzD5gZh81s4+Z2Zt3sV9CCCEGzMhW+zsFtn5namZjAPcAeCOARwE8YGaX3P2hMNsPAPg1d/9pM7sFy9I4L93E9vs8aOdfbvHBf1qQtxpXWk4XMWKJySgrIg6QqYfXE08guvFgA86qZKYivqNLowiTtqI4ON+lJn4xfo8x5m1GbTzvOMw7Bt99BoMJm0+KO9NmQwtQGpRWLRwOkKmG7tImsaLMNKkoU1exRiN8JxplDl6uihMMy24o+jC7c+djFY1dzqavREngzzIrDh7Ps3qsdfO4U68KgHfvAMYbigwcPGZ7dWe6C5n3NQAedvdHAMDM7gVwB4DYmTqA5x29fj6AP9vBfgkhhBgyKsFWcCOAT4bpRwG8lub5IQD/0sy+A8BzALzhpBWZ2V0A7gKAG0ZfvPEdFUIIIVZhKAaktwB4n7v/YzP7agC/ZGavdC+FT3e/COAiALx4cuva2uU6kkhbZYdIGifYo2pMUQA8q9LCmwvj30aZwYcX62EqqoqFZ8smx6OE1jmi8aHxQ6AxfiNrrpTDhqRCqiOJL44tPexhTqrWEyXhFQuHAyQJH5aX78H0WMqdTOkEnWVVY/oMZAyfAa+HyWTeKDunMi9J4onsW0cGxmOeG8SKaEiS+uP0vKoak1y/yeMfvgaKiMANVajqIwFn83LbouVj3xh7ZkDaRWf6KQA3hemXHP0v8g4AtwGAu/+emV0P4AKAx3awf0IIIQbH6ZmJVmEX3f4DAG42s5eZ2SGAOwFconk+AeDrAcDMXgHgegCP72DfhBBCDJFrd6ar/J0CW78zdfeZmd0N4H4AYwDvdfcHzew9AC67+yUA3w3g58zsu7A0I73dfUMW1B6wnJG55srlOP6ru8S5yCLGkqoTmSRc7WqctcdR7TM+NCt63mcsa8sOFUT5bUK6e5TPqwoulZs3FAenfZ2GbfBY0ijtXlfJhs3Fp+f0a7tr4XCAx1U2R+tNSB4t3L08rrMP4x4n0CiReTu6eRfJ+wdoLCk7fWNMYxLvCJBjtxqH3HwuZ2OvZ3OSlgtXf9FUSLuV5JrIvnWFGXRiU9Vmts4e3Znu5Jmpu9+H5XCX+L93h9cPAfiaXeyLEEKIPUBDY4QQQogNoDvT4bKpQc1RJmFZlwdgR0mnLZihXK45nmzCxcGjTEQ/5oq33OPczKTcNkdu6nLsKPvWxcGJcCxZmhsl0lw9HQMeqDh4WO+ULAYHIcShdu82hzjMEzdvVjgcyCuhxGLhs8MyTrBw93KqSBbiwHcGsa1PcXBuy2TeRMrmyMB5cjzKMAyS6Dk4o6gcRNJyUZGJr8mwHF+v/MhgFmIJuTJMuEizx0RAd4l27wMc5OYVQggh1kVuXiGEEOJccSbvTLs72trag3TaQ55lMmerJ0EMaWhCJl2yUzD8ZOqj/PRx6PaRcjsXC28tDh51Xs6FjbNRbi/nsoYDVsl4YUVORy9WF+GB/lmIQ+b0zQqHA+RQJckzTle5vdHB2ybzRqZsLe1xp5AVBy+qxhw0tmUhDTzNxyq6pBfVZ96c1ctVY+ZJaEPh9E0qQi2nw2seAYBmbMWwhbZ5V2Vn8rFhr+qZ7s+eCiGEOF+MbbW/FtoqmR3N801m9pCZPWhmv9y2zjN5ZyqEEGLPMdvKnWmXSmZmdjOA7wPwNe7+lJm9qG29Z74zrZ1vzZJWn3nL8my0Hi7BlrjvFoXTN8/2LEuglfuWSqcx97Pl3ExzhHs4dBfFvPk2u0P7FvN46aCPkn2tXJZhPWOSA2N5NnbzRrmWZd0sxIEDA7LC4VcnlL87Pq4peJhInpPJImlbI1x1UzJv3AcKmLgaMoc5pOHKQXk8Ctmbnc/RQc2ye5XVG928zZJwdi7V1wAa5+XQhoM+xcHXeOR0jb1x+m7HgNSlktm3ArjH3Z8CAHdvjbaVzCuEEGJ4GNaReS+Y2eXwd1dY80mVzG6krb8cwMvN7ENm9mEzu61td8/8nakQQog9ZXWZ9wl3v3WNLU8A3Azg9VgWZ/mgmb3K3T+dLSA6YD2yLDlLM6o9aTYvyUJ1Nm+zczCGOKQqUItElOXoLlJJq1ni4n2NZK7fKrRhxHpxdDCTPBuWZQmW3b0ejh27NQu3NWfzBmHngAf6JyEOHBhwmJQK4yCCLMc3SqCTWanjRVfsaFFe9tXX1TwcH5bZyg+2bOMB9nFZ/lIM7mLO3437mmUTAxzM0Jy/W2XzVqX2gpRLjxOmSaiHZ48TeuRwlyXY0JnazWsnvu67nkFghsV2ZN4ulcweBfARd58C+BMz+2MsO9cHmlYqmVcIIcR5oksls9/A8q4UZnYBS9n3kWylujMVQggxOBzAYgtu3o6VzO4H8A1m9hCWw/O/192fzNZ77jrTPnJGNu+4CHTIpYh5EotalApLZCFuH3GZqCAx8em3aJivjT7BC3W8a3dJuHkH2DFcbiTKwGwWjQ7d0Sh3Psd5J2SVjBLfIQ/0D9M80D8Lcagk4ESOZHfvLCk5Npk3t40Oj4/BmOTZw8xuTXJxIQG3uYLjeD+et6OUW7t3E9m7ksiD25rasnJ6/NllUu40hIW0Xb+L5Hsgk2ez0pB92NR33y7ZkszbpZKZA3jn0V8nzl1nKoQQYvi4WZVqNWTUmQohhBgk27oz3QbnujPtk2tZBTPMT34N1FmaUe7hYIam+YA8/zaXWUkmatxiTh+Hbh9JOHP3prBZtNhGKVXG0IYp5faySzi6e1kqnMTQBi7PFrN5E9kQIMcuz5uUYJvRvheO3Xl5okUHL8uj8T1y3i27m9MvhcpRnc0bPuckYzfL32X3bpXNmwReFEEMJK3P6RqJDt7MzZtdk/N5fk3Ea5/LrFnRVjRV7t5SEua2k1/vJQb4HmXznuvOVAghxDBZGpB0ZyqEEEKszvbGmW4FdaYJq8ok4x5S7iJxvWZZvXNLpNMs3KAH6zh050UYRRLakAxQJzUUC5LRDsbHC7McGj3MZrn8lsrnscQWybwxxIFzeznEIZMc4zQHQ1RhA0H2qrJoY/7uvGybzo6Px+iw+2U/oi+zQhKuLKmUaRuDRKjtSghtmNL+TBPHciWDx2xeLl8Xy6rx9cJZveHzYjk/LbOWnTtc6a7I3y3nzRy6fQJjurIPEvC2hsZsi/3ZUyGEEGKg6M5UCCHEIJHMu2NYsliMm9v6rCedN3HmMSz3RDKZs1pPkJE4trbIkOV19nFgBvpIWtWA9Zjjm7zHNMCBnYr0pots3HG5kTLggtyz8+bwhzF9WLEk2wGXVfMoz5Y7O+dybWG6cvpG6bYqz0bzFtm8zTm+E5JDx5Pj/Rt77ua9El4fTOfJvHloQ/wi5PGC0aVbSblB9uX3yOu5EkrU1ZJ42AbL5yS1F1Iu5ywH2ZfP11mSt1uHrhy/rhy6iQTMlI7dPOChaJslbT2WW+yo13CzSo4fMmeiMxVCCHH20J2pEEIIsSbqTAdMn0HN2eDoSBXoQBJOlIFZgo2DwCv3H5dki1LQmCTPeSZ5nrTXJ9M1fCGTdYFS2q1djhsKbYhBBImXbkwfXrWvSWk5z9y8oa1ygJIkHKXdw0RyrMqzJW5edrZOguR5QIEOIz84fk0n1mjS/BiAB82POas3oZB5k/CFyqGbtSXZvDxvPJaZrAtQCbbM7Q2WbuP5QZ/5jM6J+OghuQb6ZPGm8uweOHYzXKENQgghxLponKkQQgixHrZf40zPZGe6C3mD83fL7fOA7G6Ov6zkGlDKSBOWhOMAeZKb+pyOmbs2C2Jgx24hjbVkDneFs0/HwYU7AUnbQWblQAf+sRvdvRMu4xUk0RD3+P8AACAASURBVAkP9A/SLocCZCEOdWjDPGlrDnGoy7Mdr+fqhLN5u4c2xHNpQrJubGO5mO8istCGKMlWEnBw7PJ7rGTfsN4q8CLK5y3ZyVHmnZFkH6dni7KteESQhKwAgHfN1KX19MkQ7/rdtw8SsKN+XDJk9qfbF0IIIQbKmbwzXZW8AkNiGOhj8Jnznenx67biwp7cGca7BI4a7POTqevdZ5upaJ78Yi+WW7EICQDEG84R3+2Gz27Md0xs7ArL8p1HLBY+p/VE09GMxm5m407TwuFcxLoad3o8fXVRto2LcaZU5Dy09TEgMbxsRmFA4jvDcHd89SCJE0zuRJfTSWH17K6V4x6TqjGZyhLbpvPmawAojYAHSdWYNroaI+vvsx7bSMak7hI9MxVCCCHWwM30zFQIIYRYl0plGzDnrjPNpdvuElYW+cWSTVcZuJJ1k+Lc1VjN8AOOowZjtRWO5Kv2IYkB7GMqKiXpcj2rFwfnot5BYmPDT5A5uTg4H4OoJC1osG8hVfI4U282rfC402gkYvPL9VFmtuaxkgAVGWdZs5BHaZxpkFLHnp8Dh+E1y2xR5m2TfLuOM+U4wfg+uMh5ZkjiYxWPZV0cvPnz4s9uGs8zegwQz622OMFsfHX8uNjc2EeePUvsW9WYc9eZCiGE2AescskPGXWmQgghhofJgHQmsUSqrceVUnuQd2Y89myRyKocPRjXQ9JlHGeZSbl8cvap6BJlqkzWXS57/LrNpdwVXm4UZN8Dnjf8ouXjUcfFHb9mGW8cJD+uKHMYjiVLg4fkFp0FezFXJSnlSHIBkzx5dRzl2nJ/orvXeJxpkHav9nDvTiiW0KNrvGXZbJxpHEvKY0eju7d27zYXAM/G6FYOXY4FDJ9JFSmZxU0uTp4POKk4+PHrumrM8Wt+hMTfPVkk6qrjR4c47nQ5znR/ZN792VMhhBBioOjOVAghxCCRzDsg+hUH5xivZjksrjeTgJnM4cfSaRo9SM7WuGwVYFDM2HlXewUxZI7dzKWcFU5vM/LFZaf10o3LsXIUZeAJH9cgpztd2NHBe0CyYeUIjaENSZRdFXPHYQfhoFTVZ0ZJaEOME6QPy1rcvZFRj0LzUWqvCoDHOMGkyDm7d6+QfB3fcxXpmDio2f1dBHckcYJpRZnkEQ5Quvzr75puQQx9yENoNrONrWKmOEEhhBBiHRzLFKtV/tows9vM7ONm9rCZveuE9reb2eNm9gdHf9/Sts4zf2cqhBBiP9nGnamZjQHcA+CNAB4F8ICZXXL3h2jWX3X3u7uu91x3pn2kjrYC4GUbyTth2dqhG9s4y7OcniTFuaMCyRLWuCWoodhmIhFnQQxcCLmrlJs5e3k5lhhHo+ZtRBmc31NVcafI5iVn9qhZ/osVZQ5Hze8fKL8U2ElaFA4np+/VUXmJToJ9tM7tPW6bLCi0Icij7ErOvgY4gSYGPoyqgIvmOwJ27EZplyvcRMdu5d7l9zxqls+jE5qzeDm0IX62VaZukr/bJ7QhunT7ZPFW3z1J8MxZCnjYYtWY1wB42N0fAQAzuxfAHQC4M+2FZF4hhBDDwww+Gq30B+CCmV0Of3eFNd8I4JNh+tGj/zHfaGYfM7P3m9lNbbt7ru9MhRBCDJc17kyfcPdb19j0bwL4FXe/Ymb/E4BfAPC3sgXOXWealRZatQQbD7Lm9YwXzfJOWRwcjW3cbnSSRVmTZeZ5D/2ha6Zua45wD/m6K7zcOLxnUgoL2XlEdmsOvIjq4Jj0+3kIceDybAfhg64yW0lyjPLtjCzVMTDgKhU5H1OIQ5Qyr5B0GqVdloBjaAM7YgG+KI7bObRhEZ2/tH0mumt5iEMsrVZLwM1l1SondJK/G13SWRYvT1cZzNHRnYSV8KOOeRrMULZZFujQ8v3SxF44dhO2KPN+CkC803zJ0f+Ot+3+ZJj8eQA/1rZSybxCCCHOEw8AuNnMXmZmhwDuBHApzmBmXxQmbwfwR20rPXd3pkIIIfaDbdyZuvvMzO4GcD+AMYD3uvuDZvYeAJfd/RKA/9nMbsdStvkLAG9vW6860xWxPlXrEzerJ+ups3rDwG7SFAoJNClV1kbXTN32HNJm53Hm4I3vuVcsJymVUclkB6aR5BjdvXOSa6OsWcvVQRqsQjSSEl/sJE3CBaqM3yBlTki+jrm9IyonGGXe0bhsmyahDWkJtizYGWWu6qxHMEN8H7XTlyThcCFcobbomp7RMfckq3fKn3Ph5uXzI6yz5TFNvA75miycvn2+W3pk8+6b7OtmW8vmdff7ANxH/3t3eP19AL6vzzp3IvO2DZA9muebzOwhM3vQzH55F/slhBBiuCyOUpD6/p0GW78z7TJA1sxuxvJXwNe4+1Nm9qJt75cQQojhskxA2p9xs7uQebsMkP1WAPe4+1MA4O6PrbPBTbndyrzMbjm9QD4gm52/Ec7yrAaPFzIR5atmv8Z65KmumqmbOXar9URZtUeOMIdPcPhBJB4fPjZz+gyizMqqUnT3jklbj+XZDkjj44CHKPvOk5zYSVK6DShl4EmS2zthd3GQWUe0b2yFjr/sJwt2F8fpXNiKEizfLURp9+qYwyeagximI5bIg1xMbt54rFjmnTnJxeGYVI8w4nmfXJNVWcLEhZvJs1lIw3J6/Q6mWuesuX2Rm7a3xxZl3m2wiz3tMkD25QBebmYfMrMPm9ltJ63IzO66Ngj3GX98S7srhBBiCLjZSn+nwVAMSBMANwN4PZZjfj5oZq9y90/Hmdz9IoCLAPDiya3db7eEEELsFY5cfRoau+hMWwfIYnm3+hF3nwL4EzP7Yyw71wfW3Xi/Emzd541STO3M4/VGebZsy7J50xJs7BQcxfma9romddauEcQwnzXLX123X81LF1Z0qPKJXA6gLw/IiJ7DxP3j3N7CsUvSehzMz892OFM2unkriTEEM7AcWbl7o3RJsmZ0987oRIvHqjpYs+YkE2dZNZZyayndFqXdaj1F/i47dJtDG1jKnSVSbgzD4HOnzlkeNbbFc4Adu/Mkm5dDGw6S8JY+jvtiucy9mwTUrLONxQ5vwVSCraR1gCyA38DyrhRmdgFL2feRHeybEEIIsTZb/43RcYDs/QC+wcwewtJ28r0U5ySEEOJcsV8GpJ3csHcYIOsA3nn0t1VKJ1wuU22qMr0lg7VnRbhBudyYpgsXbCUThcH0o2Z5lmEpN1suylhVyTWWv7xZku4jQ+fEQfkl8cSu5DcaeB8DF9g9Og4BB1VWcSzdRtLggg7I3JtlxCgBj2n7B1ySzUIpNQptiO7eK1S6rQhtoPP+Cn0LLMKJWJVyC7Jbn0uCZe8o33I2byw7x+5dLkk3TQIvouybZfECZYgDP5aInzvnOpePPoomzGZ8jR6/5kdBxXytbt7uy+4zW8zm3QpDMSAJIYQQx1hdT3fIqDMVQggxOHRnukf0kU+yZY3Xw3JtMm/mkM1CHCyRYGv5uodjNnEVRmk3k3V5PSx/8fvquv3RiEMb4jTJ5+WStB4qgZYMvI8BDyzxTaI8zCEJ7PDumL87J+dxVToshjaQ5DlGIuWieeT9YWNLTVsebyQ+7+IvxSjtTsmxG6XdzL3L7TN6j9OiBFtzSANQunn5MUD83LPznM/ryrGbuPptfvJ8J02XbY1Na807DAyLPSpsdq47UyGEEMPltAIYVkGdqRBCiMEhmfeMkkkvnLfbJzszSkGc5ZmFNjjl1LK7d1VWzdTNHLssfy0SB3MOvcd4PDh/OJzZfKzYiRwNo1xKbTFuLr81CXbrrDwbrzcrzzZOyrMt26OsWR68Mre3lDVH4RhMuUQfZ/XGABDafpR5xxxikYSSz0dJ2ELVFnJ7E/cuAFwtJHIObWh283KIQzzvp1UZvmYpN55LU3Lv8vfCqpndTFc376bk4WrZDYVBnDXUmQohhBgkihMUQggh1sAV2nB2WFUKYadtDGoYVzm+zZJnFuLAkhpvc1VKdyLtT48gBk/k4iyPON+3crp09/LxCLJqy7GK62W5PEq7McCB52VpcEyyYizRlpdny0uwTaJ0SV80sXzbFZK2R+FNG22Djb6LRXNoQzSEzFu+5+JdBT/7KkIbRiyJN5dg61NmLUq7mXsXICnXm8+BKUmnnpzLmau/fmx0/Lr6/khGHeyfQ7cfujMVQggh1sBNBqSdU1U1GDe3rRrFVRcA77Yct/Ov1Tg9S4w6AN/tleuJJpt1fs11NRnx9udk6snuPtlotSrlLpABKakaw4zCvHSTVBqHaL/noYh1VlFmuZ4wXjWpKMN3tDw+NFaGGVeFxENFFzbqjJJLnQ9PYUBqrtTDY04zSa4aZxqOwZT2LUYmZhGBQDl+NGtjA1JmAuPrMLv7jNdIVQw8qQxTq1fHr9uKg6/KPsYQznVnKoQQQqzO8pmpOlMhhBBiLVx3pvtBH3k2na8yGjSPL6uMB4l0yiwKmSgZ40f7PTaOFwyz9ijcXRh1ElmXl2X5Kxtb2idOMI4zHfMYx4MgebaYpVL5OpqMSLodz4/3Z0wnwZirpAQZcUqf3eGo2SiTGZImJGteCfNynKDFqEE2IJE6Gw04C6eqMUVB8NyBFO8qOBauNBk1m4oywxHQvTJMbUDi2Mgo2ZNZKrYlY8GrYuCJXJuNK20bp951/OhQpds+7JObd3/2VAghhBgo5/rOVAghxDBxaGjMqRPljUVzsYyV1wn0K+4bySLFuK2SR4PMWRUZDxFf7EjN9Ic+hbujtJvJukC579m8fQqFV+P4ivdFsmZRyLxcj5FcPIvVPJKoQadxptHlyWMTeTpWlZmQlTOOeeTi4Jm713gsaTgGHDUYpd0rfD5ksnvi5jVvfnwAlGNSq3GmsdpLNZb0ePpqVQmG4gWLyjDNjl0eV1qNJY1VYzgKMhlfPZ81t/GY8vJxD8q2FaP/+kSX9mEb36H9sb3qTCXzCiGEGCSLow61718bZnabmX3czB42s3cl832jmbmZ3dq2zjN5ZyqEEGK/cQDzLQyNsWUE2D0A3gjgUQAPmNkld3+I5rsBwHcC+EiX9aozDWQRX3XB7WMyiQYoJWGWhy1x+LE8aok8WkYUNBcHb5NVV40BzBy7mZScyd5tpMXBp8evJwflcuxEHgUnbravHCU3KWLmymMeK8os22NR7yROsCVcILpbuZB4dLZeMbq0w+4d0iFn2TdG600omGEUt9ny0RXFwWnm0rFLARO9QhuCK9hJEg7HjodYcKRjlH3ZzTufN7cVrnWqGpPGCVaxoqFtjVEGnUcg7Enlly3JvK8B8LC7PwIAZnYvgDsAPETz/TCAHwXwvV1WKplXCCHE4HAYFhit9AfggpldDn93hVXfCOCTYfrRo/99DjP7CgA3uftvdd1f3ZkKIYQYJGuENjzh7q3POU/CzEYAfgLA2/ss17kzNbOXuvuf9tut4bHqoOY64zd+yOSqbKkeERkXg77JZVpJqcevuYh11NxYbuhTgDsbhN4niGEWJC+WcieZczFzllZvLOwrOW0XsTg4fXZ87EZhX2cU/hDdvVw1Jg7mn9ABmHEh8bDNGb2RGCBQuXkTdy8XEh8loQ1RLh2xCzf5vuJs1HG2HiI6eHk9UdpNQxsqN2+SqZtI5Ozm5c8gXmtVBnNsmzWf9/yYpg5oOX5tSaBD/l3TvFxb2z6GOGxJ5v0UgJvC9EuO/neNGwC8EsDv2vIcfjGAS2Z2u7tfblppH5n3X/A/zOyreiwvhBBCnDYPALjZzF5mZocA7gRw6Vqjuz/t7hfc/aXu/lIAHwaQdqRAh87UzL7JzH4EwA1m9oqjW+BrXFzlnQghhBAZ10IbNj00xt1nAO4GcD+APwLwa+7+oJm9x8xuX3V/u8i8HwJwPYBvwVJH/utm9mkAfwbgr1bd8K5YT8rtNm/m0F3OGyTYHhJwnRMbtlG5YJudrX3o6tjNZF3ev0r+St4zH59suUVSHDx+BFMwJJ2GFfP7mIyb5b9pkIQnJOtm7l7O+C2yeUmCzty9VW5tDFTg/N0IHeIFBRhMgjw7oaDnwkHdcprF512csRq/8K5U+cPH0xzaUAUzJGXWyiAG+uzoM4iyPEv0RW4vHav5NMjeifseKL8H2iThsi2fPstsK7TB3e8DcB/9790N876+yzpbO1N3/xSAXzSz/+DuHwIAM/t8AC8F8O+7bEQIIYTog8POZj3Tax3p0esnATy5lT0SQgghoBJsZ5Jahgk5uVVuL0ueaJw3k4BZJoqyLztUS4Upd1mW62w+WTPHbibrAuUxYPkrvk+WxPswKi60ZnfzjKVT2p9peC8jcvNmub3R3VvJf3NeTwxtIKdtdOiStp+5e0cUmhrDJ+Z0PGbxfCXptPq+CotW2bzFelrcvFHmTUIbZpVjN+T2kqx7lYIZroZjkEnildOX3LxR6ubPcpGc9/FcmiRBDMym8nar75MVwxiGKh3vUzavOlMhhBCDw5HXWx4a6kyFEEIMEt2Z7gnr5FOuKouwHGqJ5MlycZQrWaocF5JadzaVqZs5div5Ojqh18jmjW90zoO8Do6P1WRKzk06QuPQPGX5Ojp2x81SLsu6C5o3yojsni3colyOjCXhIGuOqG0SJFDjIBFrlmQPaTo6b9nNO24JaojEkHL+UlwkpdNiMAM7dLP8XZZyYxgGhzTwZxDdvZVrO8j7nmTzVi7+6lpHI/mjoHy6iU1Jt7ye0yvJNmzOdWcqhBBimDhMBiQhhBBiXRZ7VIvl3HWmq2ZZ9nHfjWjA/rgYrM1t4XVV8qycNyvBFi2Y4x6/5rL8Xd5GnyCG8TTM2yOIIXMyLsaZxEjLxaSGg+bjCADTIPcbybNR1pvxe05ye6vp8LmPRywzh/WQVMnnUpRrJ547XRtpOT0mQQavSpdZdPOW7yP74qvWE6XcpKzalL6isvzdzLHLIQ1cTq9rmbXs0Qef55nLPyvP1odVQ2mG6t5lWI4fMueuMxVCCDF8HHWBhCGjzlQIIcQAsaJo+9BRZ7oiUY40kmNr6TJpWzRLP+wGjFIQBxFE7Y7l2ej6zEIaeNk+QQxR1gVKyauPlFtnDp+8TgDw4hg0Hw9vkdTisZzR+4ghDkaBDllubzY9H7Mcebz9q7NSqh0fsJs3yLzk/DVvDlSI5dmuNB+q5f6EfxyQlBvXyy5cJgttiNvg9US5usrbrdy9zWXW4rwc0lCFNhQl2Di397iN3d7RKT6Zcluz7GtJNu86bt4+gQ9D51rQ/b6gzlQIIcTwcIU2CCGEEGuhO9MBEMMYsgHGWd7usv3k1ydNp/uTuVmLbXC4Ae1PUJ+qIIKD5hHhbdJusc2Ojt1M1uV5q/dVbKPzrlVEqbAWHGNoQ9v7P156QZp0PHY8YD/KfweVQ7g5x7dydIdtjki65bCBSRjdf5VKl0Wn75Rl7/C2Dulgsex7YM1uXg6D6Aqvp5BgwdJtc97ulBzLs0LKbQ7DmPUpwUZSbnEOJNdEVXIteWzTJ7f3LEm3Z5kz2ZkKIYTYf2RAEkIIIdbCJPMOiVUHNfeBJRt26pXl2prbKvcsu1cL2YhyWYPsO0/DDUpYyo30CWLIHLu1Kzgut/rFErdJamgp+ZEEXmWohmXnJPHNgoOXy7NFd++UHLqcm5vl+EZ374yc4SxHTsLOVi7yKPtmRlv6rFj2jXcD4yTpeZzk/QKleYS/FKPsy3m7qZuXPugr4fhcoWc6cd4ZHccqfzcJbYjnBLt5DxPptn5sE9t4XnQmc/6uus5s2dPK4nUotEEIIYRYG7l5hRBCiDVR0P2A2VSWZebCZbdmlHe4DFPpbCVpkOeN26xOsrhs9xMwcxxWJeF6BDFkjt0iqCL5PFge9iqoopibpo/ndZJnD3jOIGPxvma5rFlub+buzXJ7uVTahN5zDHUYTUgSDpusSrBlx44+y0ksnUbu4pjHO+9h7OXc3sLNmwQzTElj5Hlj+AKHLeTZvCz7hvUkbt4qrzp5ZMHXehF6UoU2oFNbG/uSudsFybxCCCHEurhJ5t0X1jEB9FlvWTWm+S6WjQ/VHW80yvCdR3Fn1v2WoU/h7mzsaBWdltx9ppGB2WdQGXfC9pI4wcm0bOGa8NweuRIiA2d0h5tFDWaGpCxqcEp3gpMxmadC+2jRPM6Ui4xfCau5rqWqVTQLjagYeNz+uOU8i+OAeYhD3AabjLJKMFO+48/iBBPTV13kOy5H6wl3qnzuHsyaz/NqnPaGYgBXrf4y4hO/YZ1DYXlnetp70Z2dFIszs9vM7ONm9rCZvSuZ7xvNzM3s1l3slxBCCLEJtn5namZjAPcAeCOARwE8YGaX3P0hmu8GAN8J4CPb3ichhBDDR6ENJa8B8LC7PwIAZnYvgDsAPETz/TCAHwXwvTvYp07klRy8sa1ez8mvAZZOy7Y0Bo9cNFH2rc1Jx7CpiOlqMspk3eW8YZvJOFOGt9kZlk7jeuhYsTFkFLSkynQV1ssxc7Got5MZaNVC4mxA4vGRUU2estHNmgcEHoamK/T+D0YswYZ4Q5Jyx0XVmO7ws69Myo3TV1hyZUNSNA4llWH4OHJx8GgQY5NRjBA8qCrBNEvA2bXO52C63BYk2CHKusy+GZB2IfPeCOCTYfrRo/99DjP7CgA3uftvZSsys7vM7LKZXX7GH9/8ngohhBgMi6MUpL5/bbQ9ejSzf2Bm/87M/sDM/m8zu6VtnTt5ZpphS0fDTwD47rZ53f2iu9/q7rc+275g+zsnhBDiVHAsFY1V/jLCo8c3AbgFwFtO6Cx/2d1f5e5/E8CPYdlHpexC5v0UgJvC9EuO/neNGwC8EsDv2nJ824sBXDKz29398g72b23aqjxEGbGWQ6NMROupxqSevM6j1uOXPaTSTFbNHLuZrAu0jTNt3kYfysizZtsfqYg4qH65NleN8eCAnI9ZSg5VWkgCXrWQ+IxduOTANIuSMI+VDEXOOeovHKtDUoPZ2RrHmbLsPPfV7JV8txDlu9qx2yzdchWdaTGWlCMC41hWPuZ07MLnN581u3kPq7GkUQIumk5wymfRgye/7sum1jMI3Lb1zLT10aO7/+cw/3PQYYjELjrTBwDcbGYvw7ITvRPAW681uvvTAC5cmzaz3wXwPfvSkQohhNgOfcpHEhfMLPYhF9394tHrkx49vpZXYGbfDuCdAA4B/K22DW69M3X3mZndDeB+AGMA73X3B83sPQAuu/ulbe+DEEKI/eKazLsiT7j7WkMs3f0eAPeY2VsB/ACAv5vNv5PQBne/D8B99L93N8z7+r7rTwcjU9soLRbeva2UJ8lVmUin7BaN87L0w/JkzAHgiLw+QQ2RTcUA5pVhmreZHas+ZPGK7DyesUM1Hleadx7CFtjJeTV8BtHZu1yu3JtZcPuOksLh/LFW1WeSCjMxUKFOpogrKZvGtI1FkHnHyQfUp2pM3RbDFlg+D20sASexgJWbNwltYGd2uZ5yX6O7l6/feE5kVWIAuiZ6BTo0T2fffaKRtkePzL0AfrptpaduQBJCCCEqfPlsfZW/Fj736NHMDrF89FgopGZ2c5j87wD8v20rPddxgkIIIYbLNgxIHR893m1mbwAwBfAUWiRe4Jx3pn3CFrK2rOA3wMXCuWJIkKJYtkuycauKFGHevLpKSZ/C3X2CGKK0mkm5ffJLmTh+v09WMe97zOad02OAmL1KNayL9fBnx1JhlBXHRrLmJOTUVgXIKbRhumJYaZQD6apfUPhDdDTPKQgiSrt9qsaw5BsPF8uzUdqNVXIAYEofUAxjyCTgaeXQLaejlDufUj5ymD6YZddE0ZSe29sq6n2WcHS6y1xt3S2PHt39O/uu81x3pkIIIYbLPgXdqzMVQggxONzrSlpD5kx2puVg/s2vs62NJc84zdLPvMj45ezX8mdZNDZy8MAilvhK4rTSEmfo7tjNZN3lvHGbPZyLvS6e5mMXP3c2PB6Q7Ds9aB54H13TXKptEdqujviza3b3TtmlPQvvg9rmdA7EUAfKbKgCFhqhAzKmMm9xbB9n887DNqxle/F5F8t1McShCmJIHLpZMMN02sPNy9fhLErC9BggCRkpr4mWbN4067t5G32yes+afKygeyGEEGJN9inoXp2pEEKIweGQzLu3VAEPK8spLNM0l2uLsqbR0/Y6m7dZfovZvH3KmGXvg12vUQKt5a5y2XEmacX3vIb0FLfBTtsitIFkb5Z9o5zOQRlRvub833Fom5A+O6UPwUI274gk4Zj5O2M3L2m5sZkDHaas+zbA0tmYpg+C7MtSbizB1icnhN28cR/YhVuUTksCHQDKNWY376zZzcsxdVHa5RJsxedMjzPGiUM3y9ZmVs2ozpbbR1l3n1FnKoQQYnh0C2AYDOpMhRBCDA4H4BuKG90F57oz7SODrJqduZyOEiwHOhy/Zjdgnc3b7NitZd9ubCpTNx2gXpWial5uVfqENrDs23kbY5IKwwFy0o6vULm2KN+yzDueBIcuP2pgKbdYdlWrOn9D0fsKdwPsEG7L422iDm2IgRfNjl1uu0qO3Ssh1IGdvoWblxy6PD0rpFzKTo6PV6rAj+Y2Y3d+R8dua0nHcyTf6s5UCCGEWAdfqwTbzlFnKoQQYnCsWYJt55y7znQb0m6fdWbSz5zCTikWNXXpZkENGXWZtfA6cey2DlBPHLuZ03d16P2H9c4Oed7m45w5qLNABw5X4ONTZPNyoEMoz8ZuXlzNpFxyDHPiQ5wzCVCY05s+CK5lJ1k3flx9QhuYWZFr3OzYZacvO3bLYAaSgK9GCZidviQ7h3045NCGafP5akl5tura2sp5f7Zx3ZkKIYQQq+NQNq8QQgixHm4KbRgSLKcskneclVJbRx4u5dGybZ64Xnne0t3bHPCwSMbutwU6dHXsVnnE1WD2bpJWtj8siTMe5EheTzwGk6skuZLse3A13crxqzTQgZzYpM7GYAYu13blavKBHZZvbBTCB1jVTd9GIQnn4Q7uzfm7ZV509y86DrlhqwAAGExJREFU/lJMQxuSIIYsmGHGDt3MzUuOXQ/tLNdO0tCG8BiAzrNVRwCsk8Xb9XtKMvPmOfOdqRBCiP3DITevEEIIsTaqGnMO6FpOCQDGwaU7T9yzLC+xjBfDBjhftoQk4EQuzfJC2ZEaJa7MobvcZnjNDtlkf9pKxEUyB3MWYpHJvpm0XpVgC1LumGTDCb3prFxbhAMdjAMerDnHN8q3LPlGV+SEwmc5ZSZG/FahDaMVQxvoDiM6inn7UcplWfdKUmaN83ej7MvuXZ4+COf2AbUVZQlpX8vQhrKtT9jCpkYZnCm8zkkeMupMhRBCDA7JvEIIIcS6uEqw7Q3tmbrxNUs2zXJXms3LTt8ip7Zs61NWbB6creMeJ+CmMnWz8lNZUAWzqty1YPdsLEnXklsc3xdL61mgQ3QwL0j+dHJ9zoMky+XaZsFpWwU60P4UKmdl380cu81tnDITS7DxI6soq7ZVfMtCyhcd3bws3c7pXIrBDJzbOw/LZu5doPwsWbIv25qvCT7P+RFK16CXVR26Zw2H6c5UCCGEWJd9qhrTrZqwEEIIIRrRnekGaJNhsvYym7dsq9y9SVmxPgPoi3VuKFM3c+zWGaXx9aZkHHYwH7+eHpIkTfPyYPtVtumJQxcoQx1q+T5k0dKBtDGvN/z+5czhIPtyMMRBuNIXVRu7e4930DhzOLp7+zivSS+O7mKOjIvSLsu6lWM3CWaYJm7eAw5fWDSfr+Pk8cY4Oc/zbF5+TLJabm8qF+eW/+HjCroXQggh1kJu3oFT/FqrqonQvB0jvvoYmVJTQkvB7VgAe1TdiXYb/8emIqaryaht7Og4PT7NkYWrwqarOD704Gp5bDj6b3YYzFtcoD2OJV1xDCpQmlp4DGokKyrO8JjU8qlNVgCc7kTpVjlGCHJx8o2NMw3TPDC/HB/afCdaz0tt4e4zMxwBeWRgVjVm1aLea0UGrnjHuY9Gpn0aZ6pnpkIIIYaHLx8HrPLXhpndZmYfN7OHzexdJ7S/08weMrOPmdnvmNmXtK1TnakQQojBcU3mXeUvw8zGAO4B8CYAtwB4i5ndQrN9FMCt7v5lAN4P4Mfa9vfcybyRbckeXQ1HQGlcqapVJCajObWxPNmVVaXczFS0nO4m5eZjTttMPc2SY1xtmyQe31d2zPvsC8uIxdYT6XZBbdMe40UPgpOHqyNNQiUY/rJhKTcakliCjSa5dYqDR/mOpdy4f5msCwBXrzaPJY3zHtLnwZGBmcxbVH1i41IhAbeMve48znT154T7KOU24sB8OzLvawA87O6PAICZ3QvgDgAPfW7T7h8I838YwNvaVnquO1MhhBDDZM3QhgtmdjlMX3T3i0evbwTwydD2KIDXJut6B4DfbtugOlMhhBDDwwFf/S79CXe/dd1dMLO3AbgVwOva5j2TnWkWM9d1uba2Tbn4LMYAcuFwNMtxvI0+77Np34DuMYBt7zkriN5VxmqXrLL1hEo9vF5aLrp9qzGp4ZhnVX1qWZfiBQvnb/N+z1vGq2YVZ2jvaPr4A+Ff+zE+cNl+fDLlbt7uX3Spm5ekvBgvmI0dBUpp98pVlm6P2/jz4XMwq/4SxyHX40zjOsrl1nHsiq3yKQA3hemXHP2vwMzeAOD7AbzO3a+0rfRMdqZCCCH2G8fWnpk+AOBmM3sZlp3onQDeGmcws1cD+FkAt7n7Y11Wqs5UCCHEINlGaIO7z8zsbgD3AxgDeK+7P2hm7wFw2d0vAfhxAM8F8Ou2rCH8CXe/PVvvme9MWT7pJft2HBxdy5bN0XZ14fCwFMfaHTY7Sz1xsmYkeQHLLawYA8iO3XEy7zaKJNefa/NA+ykd1yj7csDDqjGNzKTjIDQOlGAOiwLg3X+2zxfHyx1MKMRiUe5clK8rmbfFwdu4fW+WeXlgfhwnmMm63D7hAu2hLXPv8nQlCcfC6omUm0UE8vQ6VWM2F8E5cLZYHNzd7wNwH/3v3eH1G/qu88x3pkIIIfYTW/HOdLWffOuhzlQIIcTw8NXHz59Gxr86046slaXZNad23CwPM+z0LeTicXMoAdOnUHcWxJBXyMjX27RcP5qPHUvAXCUmyr7s9B2HIvBc7SWuh9uqcA6SB5uYj7OQhpJDcuxOY2hDVRmmWxsA+Cg6oTkreDO/+YuqMXT3EaXb2bRZHgaA0ZWQeUxS7uGVKN2SBFxJucev+fyYhGm+lvqc57kkjJU4yy5gQ/tjqSGhzlQIIcTwcGstzDEk1JkKIYQYJDy2e8icu850GwOn15N3mkMb6gzXjhJbD6l0VSl3HYfuqoWQc7oFOJxElPXY6TtJCofH0m3ZfH2oXcndZd8Iq8pXwmITuuo5ZSY+bmA3bzzOlSJNZGMEi9AGLjUYluMs3syxy9JtlHa5GHg9bwxtaD63+RqN7t7tPQpqXo8YDueuMxVCCDF8zOs6xUNGnakQQohBIgPSntAmn5RSJrtFo1SZu3C7zlsFIbC7NwkiyMqRNe3LSXTN1F1nYHnXUlR9yMM5usvlmdM3l3K3E/ZQE2XO8psm/oofkWN3Fhy7c24jKXcUHi/U2bzHr/t8dCz5ZqENcWwhl05jCTbKtQdXqCRdEsTAMm/83PlzjstyW3lN5I8+tnWun1XM9yug4lx3pkIIIYbLqqENp4E6UyGEEIPDPB8HPzR20pma2W0AfgrLUOGfd/cfofZ3AvgWLIMrHgfw9939P+5i37qyKcmmjwTMkmOUcuccEhDdkeSyZFdusVyPAIVNlZ1rmm93rCb7Vpm+Yd95MP9kS7Jv3OZsXn7Qs4Pm9xGHGEwo/3dGoQ3z8eZDG7I7jAkfu+DQ5aERnLGbOXajPHvdZ7u7eTl/N14jfC1t6prYjsN939mvcaYd47dXx8zGAO4B8CYAtwB4i5ndQrN9FMCt7v5lAN4P4Me2vV9CCCEGjC9/WKzydxpsvTMF8BoAD7v7I+5+FcC9AO6IM7j7B9z9maPJD2NZrFUIIYTYC3Yh894I4JNh+lEAr03mfweA3z6pwczuAnAXANww+uJN7V8jXX/htMs7nszbZ3vN7uKCXmXMVpVnuzsXt5PF20bzMedghFj2bFaVvTuGZfem7QHA9JCOz6irK3h1oiRm5NidhO2zHMxSapSBF0lowzoU+8plCYuSZ+RwJ9kvc+weBLk2K7kGlJ9tJdmHtiy3t4/DvVepwR6p7WdJIjbIgLQyZvY2ALcCeN1J7e5+EcBFAHjx5NbTqLIjhBBiF8iAVPEpADeF6Zcc/a/AzN4A4PsBvM7dr+xgv4QQQgwUVY2peQDAzWb2Miw70TsBvDXOYGavBvCzAG5z98d2sE8nsq0sza7ZvLXk2XzzzfPWma5Ny3VvX1UC7rOerus4iew9x23W8zUfV5ZgM9m3/NXc9tk1S/Sryr78RTM9aJ53EdwRLJ1NSMqNrvFFlQ+9GeK+V2ElYf/qLNzmeQ/YhVs4dJvlWd6fg6u8zbivZduq3wPMqkEmZ0nWrVBoQ4m7z8zsbgD3Yzk05r3u/qCZvQfAZXe/BODHATwXwK+bGQB8wt1v3/a+CSGEGC7ZsL6hsZNnpu5+H4D76H/vDq/fsIv9EEIIsR8sQxt0Z3quWVUKqiWbPrJv173j5bpl6PL0qlm8J013XS5r7ypzL+l+XKMcmG6jkoOzbbDE2NzGxPc8PeTWWM6PypoFh24t69Ijg9DO4SDF1loCHDInZhoWEmVeDldIZN/MscuyLgczxHbeZjy362zek9fBbTy9TgjMmZZ2iX16r+pMhRBCDA5zKAFJCCGEOE/ozrQjPHA6c6v2K9cWJSSW30DTHSXHhH4u3Gy+1ctNbSfEYfVhx3GbmXuXyVy4o/TzyWTWXC7OZNco7c65fN9Bs0OXy6ONwk9slkPz0nbNZJ85y8FRymXHci37dnPsZrIuT1dScnD3rvPoI//OaGafHK2bhgM9how6UyGEEMPDTQYkIYQQYh3MZUDaG6rMS3JHrp7N2zydz7u6ezcPMGhuy+btE7bQT+bt7vztTvfc3OzYsfwX2zIJuC3sIcv1TSVQPq4hmIF/tcdt8jozCXhCzonsXFpsyGWRhzYcv+b3WL2vmM3LmboxmzeRddvmzTJ1+2TzZm2rhjasCn/GQ+20tpWA1KEs6NcC+CcAvgzAne7+/rZ1nuvOVAghxEDxuujAJghlQd+IZeGVB8zskrs/FGb7BIC3A/ierutVZ7oiq//qbI4TrNmfcabZencxhi6722yPXuxmOqrvWkMllsPu8/LnmrWNxs2GpNGi2czGY0fj3SgXB+fKMHFZvotd0fdWMe54fvDQiGwsZzYmtc840/outtv21zHl7YLymji9/eiKYWtB958rCwoAZnatLOjnOlN3/9Ojts73xupMhRBCDI/1npleMLPLYfriUdUxoH9Z0E6oMxVCCDE4DGsNC3rC3W/d4O60os400M80kxX8zirDdN9mtg+7GGeata1jTuqzP7shG+fZbbm2McKzw2ze5u3PF83HmaXlOAY1l4CblwOAcTTjjFf+MkvJxg9mEjAbUuIzNa72Uozh7jHOdFNxoH2KgxfrSca074pVv1/2hE5lQfuizlQIIcTw8K1VjWktC7oKihMUQggxOK4ZkFb5y3D3GYBrZUH/CMCvXSsLama3A4CZ/ddm9iiAvwPgZ83swbb91Z3piqw+drNsi/ISj03MHardt9913+q2zVSCGbrMW1afaXb+LsbNUi4vx59ldIRm26j2rYqUPF5vXeS7eZ1xmmXdyQ7GlWZk0i1/MWbSKd/FHCTSbV79pXlebusq3fL0Ouf9Nq6ZdGzxaUm+WywO3qEs6ANYyr+dUWcqhBBieCgBSQghhFiPpZv3tPeiO+pMEzYVJ9hVns2i7JasXhnlePvdi0/vpm1bQdbrH6ua1YqKA93lYj4eLBdHiTiTgDMXMF/0mSSckVWwAboPuO9zfrCUmwUzZO7zWq5dfz19qsYwq1aUOdNsUebdBupMhRBCDI59uzOVm1cIIYRYE92ZdqSScMLAapbbcrlnM0WsmU05fbtLUW3r6SbPbOuXZxZSUMvw3TKP+xRr7+P8zSTgLFc4k2fzNnq/lGmbybceC4c3z1aRjRfM5OA+ebd5Me5yuT5u3iL8YUV5uG1/qgpWW2YvqsbIgCSEEEKsiTpTIYQQYj0MJgPSeWfIBcA3VQ5tW67g06F54H0mAUfq5TYT/pDJx31k53yd5XT2pbCtTOimefucZ9l110+C7ePY7dbGbKssYfG4Z5a0rXENLnbVa+jOVAghhFgPU2cqhBBCrI8601NmFxXly0zOXArM96G7xJitZ3Nu3u55vF3b1pl3F2SO3a6fT5v8uYn8X142l3Kz9XTPBq7X233ejE0FgHQua9YjUKHPerqWbjtpuml/NnV9sBybOYaHWHLNFNoghBBCrM/QfnhnKLRBCCGEWJMzf2eaSmE93n2at5u45k5atqmtXq5Hqa4Vf8Ftzt27nW1ugj6u6GzZVSVgXjbbnz7l2TYl5fbZxqZYNdRjGxIwT29OSk7m3XFIw14iA5IQQgixHnLzCiGEEBtAnemA2YYbcVVZt+/+7MLN23V7Q5Z116FP+bwsC7ePRN/HFZxJwqtsr41sG5tiF+Eg7RJsNxf7WjJvzPNe0U28Kfp8Z52W01duXiGEEGID7NMPcXWmQgghhoeemZ4dcidnnK973m7X7bVts23ZIS13HpyLbZ9dV9d2u8zbfG51d+zm0tmuZb1dOMrbSrllbZkEmy434/bNBzP0IY5eaBuBMARkQBJCCCE2wD51pgptEEIIIdbkTNyZ7qIkUL9fSNuRfXfNpuTZoTny+pRZi8egyjpd0QG5atuS1cqs9WFbQQ1d6XO+bCrjN2tL502ukXWk5Yyu86aPGnrk9p4WknmFEEKIdfFhdvJNqDMVQggxSIamamWcu840yht1eavV1tkmRWRuzVwq3L7ctq2TdZ/kmTZ36zVWlYeZPo8ldlEeLcscPm22FQ6yKel2K9vf0rWTBpBs6Zxch23KvGZ2G4CfAjAG8PPu/iPUfh2AXwTwlQCeBPDN7v6n2TplQBJCCDFIRvPV/jLMbAzgHgBvAnALgLeY2S002zsAPOXuXwrgJwH8aOu+rvIGhRBCiG1y7c50050pgNcAeNjdH3H3qwDuBXAHzXMHgF84ev1+AF9vZqlkcyZl3kyyOO1s3n7bGI7c1sZ+ybolq7py+3w+q2Ynb0Pq34X7fVusakhZ51rq+nkNIa+663nW9j04iOt5PZn3gpldDtMX3f3i0esbAXwytD0K4LW0/OfmcfeZmT0N4PMBPNG0wT2+rIQQQogTecLdb93lBs9EZ7o540UzfaqJ9FnPqvuzLQbxi3RAbOrzWv24rnb3m3J1tT0ZOrs4d/dpG7mqsvqyu2RLx/tTAG4K0y85+t9J8zxqZhMAz8fSiNSInpkKIYQYHFt8ZvoAgJvN7GVmdgjgTgCXaJ5LAP7u0ev/HsC/dvf0mcuZuDMVQghxxthSaMPRM9C7AdyP5dCY97r7g2b2HgCX3f0SgH8O4JfM7GEAf4Flh5uyk850G2N6IqchSWxDfpDEul/o8xL7ylBk3AzD9q4xd78PwH30v3eH158F8Hf6rHPrMu+2xvQIIYQ4w2xP5t0Ku3hmupUxPUIIIc42+9SZ7kLm3diYHjO7C8BdR5NX/tGT9odb2eP95wKS8VDnHB2bZnRsmtGxaeavb2Olf47fv/+HYBdWXHznn9VeGZCOBt1eBAAzu7zrcUT7go5NMzo2zejYNKNj0wyFI2wMd79tG+vdFruQefuM6UHXMT1CCCHEUNhFZ7qVMT1CCCHEUNi6zLutMT04knvFiejYNKNj04yOTTM6Ns3o2AAw3QAKIYQQ66E4QSGEEGJN1JkKIYQQazL4ztTMbjOzj5vZw2b2rhParzOzXz1q/4iZvXT3e3k6dDg27zSzh8zsY2b2O2b2Jaexn6dB27EJ832jmbmZnZthD12OjZl909G586CZ/fKu9/G06HBNfbGZfcDMPnp0Xb35NPZz15jZe83sMbOTx/bbkn96dNw+ZmZfset9PHXcfbB/WBqW/gOAvwbgEMC/BXALzfNtAH7m6PWdAH71tPd7QMfm6wA8++j1P9Sxqea7AcAHAXwYwK2nvd9DOTYAbgbwUQCfdzT9otPe7wEdm4sA/uHR61sA/Olp7/eOjs3XAvgKAH/Y0P5mAL+NZaTuVwH4yGnv867/hn5nqijCZlqPjbt/wN2fOZr8MJZjfM8DXc4bAPhhLHOgP7vLnTtluhybbwVwj7s/BQDu/tiO9/G06HJsHMDzjl4/H8Cf7XD/Tg13/yCWIy2auAPAL/qSDwN4gZl90W72bhgMvTM9KYrwxqZ53H0G4FoU4Vmny7GJvAPLX47ngdZjcyRD3eTuv7XLHRsAXc6blwN4uZl9yMw+fFT16TzQ5dj8EIC3mdmjWFYd+Y7d7Nrg6ft9dObYqzhBsRpm9jYAtwJ43WnvyxAwsxGAnwDw9lPelaEywVLqfT2WasYHzexV7v7pU92rYfAWAO9z939sZl+N5fj4V7r74rR3TJwuQ78zVRRhM12ODczsDQC+H8Dt7n5lR/t22rQdmxsAvBLA75rZn2L5jOfSOTEhdTlvHgVwyd2n7v4nAP4Yy871rNPl2LwDwK8BgLv/HoDrsQzBP+90+j46ywy9M1UUYTOtx8bMXg3gZ7HsSM/Lcy+g5di4+9PufsHdX+ruL8XyefLt7r6VwO6B0eWa+g0s70phZhewlH0f2eVOnhJdjs0nAHw9AJjZK7DsTB/f6V4Ok0sA/scjV+9XAXja3f/8tHdqlwxa5vXtRRHuPR2PzY8DeC6AXz/yZH3C3W8/tZ3eER2Pzbmk47G5H8A3mNlDAOYAvtfdz7za0/HYfDeAnzOz78LSjPT28/Dj3cx+BcsfWBeOnhf/IIADAHD3n8Hy+fGbATwM4BkAf+909vT0UJygEEIIsSZDl3mFEEKIwaPOVAghhFgTdaZCCCHEmqgzFUIIIdZEnakQQgixJupMhRBCiDVRZyqEEEKsiTpTIXbMUT3MNx69/l/N7J+d9j4JIdZj0AlIQpxRfhDAe8zsRQBeDeDMp1IJcdZRApIQp4CZ/Rssox5f7+6fOe39EUKsh2ReIXaMmb0KwBcBuKqOVIizgTpTIXaImX0RgP8LwB0A/vIcFd4W4kyjzlSIHWFmzwbwLwB8t7v/EYAfxvL5qRBiz9EzUyGEEGJNdGcqhBBCrIk6UyGEEGJN1JkKIYQQa6LOVAghhFgTdaZCCCHEmqgzFUIIIdZEnakQQgixJv8/7RLuJS17B7UAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 504x360 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFgCAYAAADgoJN2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO29e7B0aVnluZ7MPKeKm6CWiFFgw4xFaw1og9Wg4YRgCz0FE1384bQWhNNNi9Z0K44jaoyGBhI4F5QQG3tq1E+GQY1RRKLD+AzLqe7oxiaCEKbKYHSsorGrS4VCw7pwERv5vrw880eer87zrvfsZ1/ycvY5Z/0iTkTu8+5b7tw739xrr3c95u4QQgghxHAmp70DQgghxFlHnakQQgixIepMhRBCiA1RZyqEEEJsiDpTIYQQYkPUmQohhBAbos5UCCGE2BB1pkIIIcSGqDMVYguY2ZeY2b82s0+Z2TvN7H81s/9hy9v4f8zsv9jmOoUQ28GUgCTE5pjZ2wBc7+7fbWZfAuD/BfAV7v43G673IQD/wN0/bGbfCuDb3P1btrDLQogtojtTIbbDywD8xtHr1wK4awsd6Q0AvhTA/Uf/ugzgm8zsGZusVwixfdSZCrEBZnZoZp8B8HwAv2Vm/x+AVwD4d2GenzKz3wzTbzWzf2Nmh8l6vwLAx7G+Rh8zs8cALAD8PoD/ajfvRggxlNlp74AQZxl3v2pmXw/gfe7+pQBgZo8A+GiY7ScBPGhmLwDwYgC3Avgv3f1qst4HzOwHAXyju3/btf+b2UcAfM0O3ooQYgPUmQqxOX8HwB+E6acB+Oy1CXd/zMx+BsAvAXgq1h3pZzqs92uwfvYa+SyAL9tsd4UQ20YyrxCbw53ppwA8heb5MNZS8I+4+8cHrhdH6/30kJ0UQuwOdaZCbM7XoOz0/hDAc69NmNnzAfwc1nem39FlhWY2AfA81HemX4W6gxVCnDLqTIXYHO5M7wLwEgAwsxsB/BaAfwrguwE838xeGhc2s3eZ2btonU84+puE+a4H8LUA/vV2d18IsSnqTIXYgKNhKl8I4N+Hf/8ygFea2VOx7ljf5u6X3f1zAN4K4H+m1TwLwAfiP9z9PwH4eQD3H401BYB/AOB33f3Pt/9OhBCboNAGIXaAmf0vAB5293/eMt8h1ne1X+3u85Z5PwTgde7+R9vbUyHENlBnKoQQQmzIzmXeo5zSh83sxF/TtuZnzewBM/tDM3vhrvdJCCGE2Cb7eGb6LqwHqTfxCgA3Hf3dgbXrUQghhNgJZnarmX306Cbuh09o/3Ize5+ZffjoJu+VbevceWfq7u8H8MlkllcB+GVf80EATzMzDUoXQgixdcxsCuBOrG/kbgbwajO7mWb7MQDvcfcXALgdwP/ett4xJCDdiHUG6TUeOvrfX/CMZnYH1nevOMCTvvaLpl951FDOtwo/EVZTagvveDkrnxcvD8p5Y/uCjtTyYHW8X9NyPdNkejop2yZhemrUxtM4njZqs4b5GEvatonzh9KRbP+yNfJycZqXmzgf11VzW5ierlZl2yq0ObUty+lpmDZqQ5xeJG1t88bpJR3HuNyK2qrp8HoTT0U88EafwsROfg0As3ABT5M2AJiFC/wgayu/CFa0nkWYdz4t512E6YWVy80ttNG9ydxpPeGLabGk7S+P3+diUb5nX5TzzkL7dFE0YTpP2sJyk2XZZqvmafqqqab/cvH7j7r7l2DL3Grmjw5c9veBu929SRF9EYAH3P1BADCzd2N9U3d/mMcBfMHR66cCaHXQj6Ez7Yy7XwJwCQCeMbvF/9un3Qug7jCvPuH49eeeVn7yf/1Fx9N/9SXlGfTpL+Xp4zPu0zeUZ+Ynn35svJx9UWnCfNpTafoLjiNYn/rEMo71Sdcdz/uUw7LtCbNym0+cHs97SFfDdXY876GVbbFznYGuGiLriFc9OsjlwM50GrbP+8LT0/BeZqDj4cfTB14exyesaHp55fHXT5qXn8ETw/STrny+aHvK3xxPP/lvyrYv+Ozniuknf/a4gMzhZ8o2fCpMf/I/lW2fpsIznwrTn6S2T4b1fOZK2faZsH+fpUjgvyET8ZVwLK/Sty537hmxA7uOLtInhF+uT6a8/6ddf/z6i55Qtn3RE8vpL3nS8eunP7ls+9IQQvX0Lyia/vqGMqDq0S8+bv/LL3xq2fbk43kfub5c7i8Pjqcf9nL7j8zL6Uc/f/xe/vKvyvfx2KevO17useuKttWj5fH5woePj93THimP69MePp5+6iNlJ/zkTx5PP/Ez5fV5/V8Xkzj83HH7jE4X7ojf9rD9GXbAowDuHbisAV9pZnHxS0f9B3DyDdyLaRVvAvCvzOx7ATwJ66pQKWPoTD+B9Ti7azzz6H9CCCEuMtOBTyKXq0fd/ZYNtvxqAO9y958+KmTxK2b2PHdv/DU5htCGywD+0ZGr9+sAfMbdK4lXCCHEBcKwlvmH/OV0uYF7HYD3AIC7/x6A6wHckK1053emZvZrAF4K4IajJJcfB3AAAO7+81gnxLwSwAMAPgfgn+x6n5qonhssT54PACar4w9stSo/PJ72OK/TvN7cxs8do8zqvJ7wLGfJP57CrG3y61B5ts8z0vg+WLpdhDZ+DnpQSdST8KpsK45VchwBYBmfe01K2Wwefhkvs2dp9Auan7vNg+R5OCPJ86D52V46fUhth+FyPqQHZnHeJ7ABoId0u+xxflyX7WsiAce2w1lzG9D92B3wM1KaDp/fyvj8iG3lctl5tuRrNHwP8CW6Sr5P+HtouoyvrbGtDyzd8vTpYMPvTJG+gXsA3GRmz8G6E70dwGtono8B+GYA7zKzr8K6M30kW+nOO1N3f3VLuwP4nl3vhxBCiDOEAZgN+3Gf4e4LM3s9gLsBTAG8093vM7M3A7jX3S8D+AEAv2hm34+1Gem13pJwNIZnpkIIIcTecPe7sFZF4//eGF7fD+Ab+qxTnWlHbDX8F1IcebBM1tPmlo3SLktRq6StUEtJpsrcu33I9p3lrqKNlovDg1jgmdN/ouy7BMmqxXq4rdTYDsPxmk/KbVwXjtdi0iwNVsMpSEb0uGwqRyaSJ0/3aYvuWZZ1r9K8GX1k3q5ybdq2wfGIEiF/diS1x8+P5fx4PfHQmFU4J1d0fmayLz/Sid8LTsd4lnxnZHJsJt32kXFPTfI1bCDz7h91pkIIIcZJu5loNKgzFUIIMT5sEwPS/rnQnekm8kV02K1IlmkeiVTDrtxI5fzNpNTQNqvcq3GDtODAH36ZdLve5sAVR0mc1nFA8uy8cPOWbyxun/eFZfBFcPMeWOmCLZ2c5XJRGlyRjLicNcu+LDHOkqSeVPYdLAHTZf/EHlJ/JvPyF1+c5v154kFzWx+Zt6Obl4/5kuXaSfPnHB3efE4uo8ybOPV5um0EQMYk+X6Z9JHhx45kXiGEEGILSOYVQgghNkAy77jJpF0e8BwHRA8dDA2UTr1M+mHptBr0ncw7ibJm5dg91oUqhWigmbfVeTxQ5l0moQ1zJzdvkH05ZDw6eOeUzD0jd+8irCcGOAClxMeBDlEOXE7Z6csyb9gfkiNnQ0MbtiUBs7s3Bubzl9kyCdPnu4i47BMTxy6HSBThEz1kb25Ljjm7rePnx67tpcVrq/nRy4LOz+p6DkH33BYfDfHHcR1NRym3j+S7LVfuXt29ujMVQgghNsBQVwkaMWdnT4UQQoiRojvTgUyDdDvnkpN9MjmDUlbl7XJJyugGTPJmWUqOQQ11hm13MumWB6wXbS3O38ikCFvI9zW+Z3ZZpm7eypGZDMqPgQ4k3S6LbN5mWZeXXXH9zpgbu4mbN8qlV6k4byyl1ibzNi0HAIseslu8qxgazLAlNy8fc3Zbr7L83Sjn02OA9DxLysY6uXeXidO3j0O3qwRbBTos8vZTwUwyrxBCCLExMiCNhzxuK//VM/TXWTZmLLtLq6rNJFVkqvGpYbKq/JLGCXa/Nx1699mnEk38hT6l5SZkJCrfV9m0QHOc4ILe8yy5K4nxgofUFo0q1V0rVx4Jd0LVXVGYnoy9akyfO4XMSLSLOMHkjr865mk0JBuQwh0unWjR+FaZipLIQL5rLfaNiz5xFZnQbsl32CjuLjdB40yFEEKIDZHMK4QQQmwB3ZmeTTLZl8dzxem2KLAsMjBKP9m4Up6uTDShbcJtcbKKE9yNdJu9565UHwft+8Saj10xNpD2dU7v+TBIdWxAKiS+JGauKiidGJLYnBSnDzcZZ5pJwEXVGDqQT+gTJ5h8uWXjTIeajDY5HkHKrY45S7mTbgYkHs+cGZB43Gm8Jrh6VPwOaatQlY1/71NFpiunWzVGd6ZCCCHEBpytBKSzs6dCCCHESNGd6RZg91023ebYLZbrUVx4lciak1iJhTbH8ug0aKmZlNsm4w4tFh6Lg6NyLNOB9SDd8fsK72NB+vCUfkMWjsykwgzHCS6De7SSdavpZudvjLo7nPWQNbdVOJzdvMskTpDHnUZ6ybyJ0/cJHV3AQCrzFhGOfMz58wkyL1eUiVWFKpd2OEcr926vqjFoZGiUab/C4SOUUyXzCiGEEBuioTFCCCHEFtCd6Xjp40yLjt26AkMigSZuvFT6aZGJsrZZdAqy5hkmJy3GzUXHgIVMgm6bN11v4UrmnaVfqUH2nVDb0mMlGI4PpGCGIOtVwQxB4ptzNZGiokxzCABQyopVxZI+xcG7unn7xBD2iROMX259qsb0kqR7VI1J3vM8iROsirnHuMekMkwWTZm5d4HSwVsFsoS2KX/XVPGCJ79eT5+jEAeVYBNCCCG2gDpTIYQQYgNkQNoTBqym7bNl1BJJj833GGRdVo3pLv1U64lG10RusjSbN91EylApd2iAQ3VYK4n6+FfrlBqjI3NJbauqGk2Q6rg4uC3DfCwBx4LSuVs0yoqcE5sFOswyWfN6qgzz+ZDH+wRqmw+sGsP0cX0OdewOLZZOxy7LQ+Zi7plre16cS7SeMJ1ViwLYzYvGeauAGP6eavm+2QajkIQl8wohhBBb4AzdmZ6dbl8IIYQYKbozTYjljTgvIIPlHc7hbFyupYSTZ6ENSTZvsfkeMazV/vWQcrN5s8zfSRLa4NXOH38onJka17OkiuyczTsLUt0Kpb4VndFLduwGObAq21UVEg+ScBIgMCcZc5aFOAx1+nIJtj5k2bzMNgqAb+Dmja7pOpuXQxuagxmyQJR4LvOjl6okW/IdEpX2PjJuPcrg5Nf1cp03cXrscJypmd0K4O0ApgDe4e5vofafAfBNR5NPBPB0d39atk51pkIIIcbHjkqwmdkUwJ0AXg7gIQD3mNlld7//2jzu/v1h/u8F8IK29UrmFUIIMU6mk2F/OS8C8IC7P+juVwG8G8CrkvlfDeDX2lZ67u9MdyVnxLxMHmTN5sgo+65IqSxDG5qXA8qB51XJsSjzkkwVgxo2MQIWMnOLLTi+r8wFzMR5i5xeoHIiWzKYfh4G0E9I1vXEzVuVZ4uOXXL6FmXekrJdQCkr1iXYmp2+vaTc68PlPKcTP7p7+UTj8IUsmze2bSu0gV3Jfdy8iRM6Oqg5RIM/nyjT159zfAzAASDdHr3wdJ3NGx4p0UfHowVie5/HT9sKdNibRLzZ0JgbzOzeMH3J3S8dvb4RwMdD20MAXnziLpj9LQDPAfBv2zZ47jtTIYQQZ5GNhsY86u63bGEnbgfwXndv/QmhzlQIIcT42F1owycAPCtMP/PofydxO4Dv6bLSC9eZDnW4sZRbLNdnzDuXYFs1t1WDwAtJqVxvdBLOV+WvuYOwg23ZvE3bq9p65Aizy7ErrCLWNP9qjbm+LMHOKUN1FtqnNCg/Sn4L0t/K0lzs9G3O5uVAh2UhVdJ6OMThuiCJXlkUbYOdvn2yefuENgwtwXbYQ+YNX7bshI7HmbN4q1CN8Pl5ks07rzKgg0Tf4saPUm71PRBHDqCES7Cl30WFm3c7HdGpOX935+a9B8BNZvYcrDvR2wG8ptq82VcC+EIAv9dlpReuMxVCCHEW2E0CkrsvzOz1AO7GemjMO939PjN7M4B73f3y0ay3A3i3u3e6BVFnKoQQ4kLh7ncBuIv+90aaflOfdZ7LzrSrLNE2X1bhPhtYzU69KOVmeZ3VehIpdZaFNtAqV4kLd8KO2Y77VklayWD2DF5P6eDNQxtmiCXYyJ0Z5DijH5ZTdtqGeQ+TzGMubRdDHDjPtZ4OciD92s6cvlWIw0FYNpNAWTqdT09+Dewumze+T84K3kb4BFC4ebnMWpTMK2mdpwvXNjl/w+eeBTOwdJtdv/w9EJlyybXk48m+ozLOTGgDf5mNmHPZmQohhDgHKOheCCGE2IDduXl3gjrTjqQu4JYkhK6hBVkWL0/3cdOyPNmVNik30md/2InctJ42hYfzeCMW5GI+jksalL9KMn5jiMMhO3bDeuYtbtEoK1a5sEGqzJy+ALAK05Ohub3Xb5LNm3wofBcxTSTprARbHzdvmK7KrMUSbC2hDUUJNjo/4jnATvBVEtowp2PlRSnGoqlQ2q9rUd27jkjYpMTkOGTg3RiQdoU6UyGEEONDd6ZCCCHEFpjozvRMkpUzqudtbvNECssGcrdme8Z5JzxvnGLHbubKpf0bWDqtCpiIldR6ZPMW66hKrnWXsidBjpvTeip3b3TzVk7fGNrA8vCisS0LcWApN06zVHn1sLxED68ctx9mEuh1dGkvC0s5tQ3M5mW9nnX5LLRhS27eRZLNu8ik9aScXiUBJxnQWV51da2Hw1VdL+Harr+Hmt29Zz6YIWNHVWN2xdnp9oUQQoiRojtTIYQQ42N3cYI74cJ1plEWmbSHvx7Pmzh2uWTSIivBVils2aDvZF6SiWaF3FQumAcodJdRikHnLA/3CJ8o2liejfvObmIOmJg0TmDqMZu3xSUdnZy0nsNE4osuz+tZHk5CHFhiLLbPoQ0sCQeZ85Cl3HkIRuASbGk2b3Id8HqGhjbwNqPsyxJwfF9cni0JteBjlbmk+fOJbmwuw7cMUj8/Blkk2byVw734HmCZ9/h1VXJtYJm1zM07oVjn0XKGZN4L15kKIYQ4A5jJgCSEEEJsjO5Mx8OuXGosxWSskmxeT9y81XpCe+0UDPIOLTeP7r8Nzs3MoVs5jzPnbceAiSnLuqxGRvmL3nSU3ybk3mV378yDHIhSU4s5rSz/HWZ5riT7LsOBd3aSFq5TLsFGkmPiXp0NLsHWoy5fny+3uB1+9pU5j2NbFkyBvHzdMnwG/HlUJdjCCVTlTIfphdP2k2uyemyTPO7p831iiZR7rjhjz0z3sqdmdquZfdTMHjCzHz6h/cvN7H1m9mEz+0Mze+U+9ksIIcSImdiwv1Ng53emZjYFcCeAlwN4CMA9ZnbZ3e8Ps/0YgPe4+8+Z2c1Yl8Z59q73jakf2B9/KGwyinCVhyt97lrj3V4VMUYmo0n8ZUu/euMJlPh06rGbzfQxFfGdaFosvOO40wXNZ3SnehDuInnsaFEph/eV7xozk1E4egdp4fBmQwsALBMDUrwbvXJQGm4O56VTpDDVDK0ow6aeDL4TjXexXG2mihMMyw4dS5qMK11PBwMQxwnG8bt0J8oFwIsKQPRZRlMaj2eOd7HzJakTSRyo0/UbC363FQOfbikycPSYnak7033IvC8C8IC7PwgAZvZuAK8CEDtTB/AFR6+fCuDP97BfQgghxoxKsBXcCODjYfohAC+med4E4F+Z2fcCeBKAl520IjO7A8AdAPCUyZdvfUeFEEKIIYzFgPRqAO9y9582s68H8Ctm9jz3Uvh090sALgHAMw5u6eGaOJk6iqv7KguppRpXSmsNkg5Lt8VyLVVjVsk402I6eVt9RJM+40OrCjMdCyH3gV3yUYJFNRYvFA4ncxAbkpZxrCBJfI5leN0sbVdjE1kqjGYlHuMYJMjDpHB4NW9WUeY6knLj84RqAHMykJHHmU7CsizPMpnMG6XmVOalajxZnCCPJZ0m0npSALw6t2NxcK4MszpeT3X9Jt8LfE0UEYHctqUC4GVFmXwbWfWZVcvHvjXOmAFpH53pJwA8K0w/8+h/kdcBuBUA3P33zOx6ADcAeHgP+yeEEGJ0nJ6ZaAj76PbvAXCTmT3HzA4B3A7gMs3zMQDfDABm9lUArgfwyB72TQghxBi5dmc65O8U2PmdqbsvzOz1AO4GMAXwTne/z8zeDOBed78M4AcA/KKZfT/WouRr3X1jGbcvLGfEGC921PVhlcSI9ZF9S6WuWVatVKE4a4+jmkm5bY7c+D7bYtaamPA4U5LNonLHDt1J5uZlqS7IvIfcFn5vsjSYFw5vrkoyp4u9a+Hw9fTkxNc87+EBnQXzbJxpj7y6aY8TaJLIvB3dvKuk4DdQyt511Zhmab2ajoXeqzjBpDh4n7HXxeOeoqmQdivJNZF9+1S6OpOcoTvTvTwzdfe7sB7uEv/3xvD6fgDfsI99EUIIcQbQ0BghhBBiC+jOdLxsSwbJCvRyNFiUdGqH3zAJOCsyzq7X4i33ODf7OHT7SLldQxuqdSahDVztJU4auVfZ3TsLH8qCZLxZjBOkbcQQh1Ul8zaHOCwrybFb4XAgjxOMTtfZomybHEQraY/LflvFwfvIvFnB7yReMTt2VVAHxwuGzy97vJEWB0+uSSCvGhPpEy2YceYlX7l5hRBCiE2Rm1cIIYS4UJzLO9Ou8kbbfHkGZsjtTXI22/AkiCErLsw/2IplE9dr26FJ83h7OHQz+WsodTWcmJtbvuko481a9jVKu+zWPCxyWamod6wmQm1ZiEPm9OVwgSXJXFcPji/Zg6tlbu/sMLRR2ELh7uVR9yzXRqrQhh6fZVeZlwImVqGKTBbSwNN8rKJLuioOTjJ8lIGXlMEcP/dl5eaN53nRVF8j4frxHi5czgXvWgCc2dojrn3Jx4ad1TM1s1sBvB3rESbvcPe3nDDPt2KdzucA/sDdX5Ot81x2pkIIIc4BO6hn2qX4ipndBOBHAHyDu3/KzJ7etl51pkIIIcaH2a7uTLsUX/kuAHe6+6cAwN1b0/jOfWc6WbTP8/i8O5IvUnm0cPrm2Z6ldNpccqzegeOXbedm5jJcJtJtNmA9UxEzx/J0kr/HaOyMki8ATMKy7NBld+80FpHmEmwepWSSh5PC4ddZ87x1Nm8If0jyZdfzBunysLx8D5bHJzDLo5PVsZQ642LgfYqDb0vmjQXBKX+3kG5J1r1yUL7nIpiBj110ULPsznJ6kH35HIiFw7PAjwWVYFsum69nDm04CNdBWwk265Gx0cSZcfoONyDdYGb3hulLR9nuQLfiK88FADP7ANZS8Jvc/f/ONnjuO1MhhBBnEMMmMu+j7n7LBlufAbgJwEuxzpN/v5k9390/nS0ghBBCjI/dyLxdiq88BOBD7j4H8Cdm9sdYd673NK30wnWmfeSNwjVH0kqZ20vLsTwb1pPJqCyHZpmyLI/OJlGqbIalp4w+wQtZjjDva9fQhgXLW5zVG6Cx/JgHyY2XY3dvUYKNJOGyPBs5dhElYHKAViXYgjzLEmPi5mUX6mGS4xvnnZEEPAkfCAc6pLDMVn6wZRsPsI/L8pdikG85fzdK1OzezbJ5K8ducVxzN290Y1fO7OL86O5a5/O+eEySlmBDZ7JRBnWJye7rGQVmWO1mnOnjxVew7kRvB8BO3d/EujTo/2lmN2At+z6YrVTjTIUQQlwY3H0B4FrxlY8AeM+14itmdtvRbHcDeMzM7gfwPgA/5O6PZeu9cHemQgghxo8DWO1onGmH4isO4A1Hf5240J1pJpGsGVYFjksmxbWskmxeloWybE9W1KLExKffUPNfH4duHyl3aIgDv6+oMrIaZFESJ5flgmTfaXBGV1JdmF44S7nHJ5Cz0zcJcagl4KxUGO17UnKsyO1dlm2z2fHBukoS8CESFnSRROfvrEUujuYRnrejlFu7d2neabNEHoMYWFqvHbvNARxlsAqtp8f1W+ZwF02pPMvfU/xYqSv9HnEN28a22ZHMuxMudGcqhBBinLhZlWo1ZtSZCiGEGCW6Mz2HGEkv08Lpm3/gbHos24a5a2sJNkqXJBN13kLz9tbTx6/bZN3MeZxtIzIhOXZJ7+tgenxgOUI2vutqPYnrkqXCWXRyJqW5OG83C3HIsnmrkmtViENzubYol84XLPMeHyCjbN4Jafbpl8JkYMBDkrF7lbJ5s9CGKps3HAM+VvH4VCXxks8rC+7gcyeGhbQGmYTvAc7zzsqusbu3azbvPnJ6d4oBvqNnprtAnakQQojRsTYg6c5UCCGEGM7uxpnuBHWmCUNLuXHJpKL0Uurwa5YRASAYMvMB4pUUN+yEzBy7baENcd40/CGTwGm/SZkrghmi5Lve5vHOLym3d0lZtKuQwMEZv4fhYq4CHWI2L4npWYhDJTkGKeuwkoCbQxwOk/JsM3LhXgnZvCzrXkEzE/oyK5atLKnkdI1BItR2JWTzLihxIzp4WdZlabuYt8rbPZ52ykqusnrD58Vl+KLsy9dEvA74XObD40X+bvM1sY3s3fPALofG7IKzs6dCCCHESNGdqRBCiFEimXfPsMwazYr1AOhmN2Im61b5u9m8lZQbXzeXZWLptJKEE4fsxGIb7VAPB2bmvO2TQ7rK5K+hoQ2kscULbb5sFlkmpJtVjswg601pZ6PsO2UJNoQ4LEnbX5LoE4MA2ElaSMBVaEMSUlC5V0MwA4UdWDgJr3oe2hBl3wOySZefQR7aUHw+5OYtAiay8AmSsnncYZR2s2NXOahJao9SLgdwFKEN9Bggnstcgi0LbWDi90tVci0JcWgLeCjaknKUfZZb7anXcLNKjh8z56IzFUIIcf7QnakQQgixIepMzwhtbt2ubt5swDXDEmxZlqlsq0qyJQ7ZLJs3C41gumbqZrIuUEq71b4mxysejymFLfBg+kkMqiDFMW5zzlIYy6xB+q9cnkEi52OzCPJw5QAl92iUdivHbnQFk3vx6qy8RGer45Nyym3L4zaWTmchq3dOJ9Zk1v0xAEvtGfGLMAtfSMMnqtCKZtmb543HMpN1ASrBRhnM8XOuAh3CuZWFgaynj19XQQxh3tPI4h0jrtAGIYQQYlM0zlQIIYTYDDtb40zPZXGs0ckAACAASURBVGe6C3mjysDskaXZVQbOSq4BpYw0o23E0AAOMOhzOnaVctscukVpqjS0obltUYU2sO59/JLzd+dJNm9WZo1lvLIEW/NgfnaAZiEOLDkukkCHLMe3zu2NkiedA8Hd2ye0ge8M4rKZu5qnq9CGsD/zQy6zFrJwk7xdoHTzskQfj2WWxQuUMm+awUxu3jK0IX+cEaf5+2NoebSh33VnQQJ21I9LxszZ6faFEEKIkXIu70yH0l4s/BgrxoXly8VxY9mv1bbiwp4YkOJdwpKj08Iv7Snd3fWp6NLHVLTMfrEPq7le/ZovhvXRGL8ZYkWZss3oTjUakvjOYxbGqFbRft5sWsnGnWaFww+5GPik+c6MxzVemTXffU6Cu21yULal40yX2TjTnGKcKRuHirGkze+Rx8tWxdLjWFIek5vdtaLZkFTFRkYDUqLWVFGDHCeYVHuxHnet2fdSdteaLVetJxmTuk/0zFQIIYTYADfTM1MhhBBiU1hlGzMXrjMd+vA+lVb6jONcsjwaXzfHmAGljFSN1Qw/4CxT4npE+2Umo7YKN6UkXa6nT0H0YjmanibJdpOwP3w8soo3LAnPwi/jlXMsYZR5myVgoDQSVTJ8lJmTijJAOUaVZc1ZMB0t6NnDLMrDThKwN58wmxiQItk406oyTDQnsQRMY2tjOx+reCxXVZWY5s+LP7t4bmdjSTPD4Lr9+PUBzRs/Lq461VXWPW+ctaoxF64zFUIIcRawqnTemFFnKoQQYnyYDEijYuj4rbZl8+K+LI0dv+ZovyiHthUXXmZjN5Mxl8U6ehQK7zN2NHPsZi7lLL6QYRduTpSHqBIMSX5xvdNJs1Q3p/cRC4ezNHhIuvMiiNTsJC3HmZILmMedTuJYUpKkgxOZ5dEoAVsi6wLl2L4ZuXk9usbTteTjTK8mUm5092buXSCXyKO0zg7qajxx+Ez4nCylXI6iDG30iGCxyMaZopHsuwVoc+ymi3Za51hYjzM9OzLv2dlTIYQQYqSc+ztTIYQQZxPJvOcQS2SZOqSB2oO8s+hRMDgLeGBZNSpMkx7Sads2i20kQQxcmaWrlMtSdsaE4+qyzyTMypFkVRhFUvw5SnfVwP8g7R5wnGBSRaZP4XCOE4yO1dmkvHwPgnY4ozjBKwcDkzKISY9C8/E9V47d6MJNHLp1nGASr8ihGtZ8XDnucZXFRhbBDEloQ/IIh9urot7JuZx992RkEvAYZd0Ks53FCZrZrQDejvU4gHe4+1uo/bUA3grgE0f/+t/c/R3ZOtWZCiGEGB2O+gfSNjCzKYA7AbwcwEMA7jGzy+5+P8366+7++q7rVWcqhBBilOzozvRFAB5w9wcBwMzeDeBVALgz7cWF7kzrwdA9iiTHQdYtMgwPwm6irbjwLMkBjVIyS1hcZLvrPgx16AKldFoVPR8Y2sDryX60zsNr4+AB+tyjBMjyaFZRZh4kYc7UZWm7yPFNKphULuCswgxpg0U2L1nDLRy8KUmlE28+kBwwMfUY2lBugx27cdmqLezDopJuJye+BoAlZxUnBcCjE5od1JXUHo4BX4fxs+PrLp7L2WMZoHT5V47/JLe3zv4+ebmTps8yG1aNucHM7g3Tl9z90tHrGwF8PLQ9BODFJ6zjW8zsGwH8MYDvd/ePnzDP41zozlQIIcRIMYMPl3kfdfdbNtj6bwH4NXe/Ymb/HYBfAvD3sgU0NEYIIcQoWR2ZkPr+tfAJAM8K08/EsdEIAODuj7n7tQJK7wDwtW0rPZd3pplTLW9jySQM5u8T/pA6dMvpGOLQJhOtivAHzuaN8luz07eNrpm6rTnCyftKy75lRdd7OEmL7bHbOnH3clmzWdhmVp5tvmqWboFSvl1Q3MEyaavWE36ps2M4hjZkub3zltCG+L44tGEVlyWnLROj4LjMWpR2s/zdK+T0rYqlF4EXWcHvPDt5Xsi8zZI9nx9Fbm91nheTRRgDf5/ER0V9sr7PMzssDn4PgJvM7DlYd6K3A3hNnMHMvszd/+Jo8jYAH2lb6bnsTIUQQoiTcPeFmb0ewN1YD415p7vfZ2ZvBnCvu18G8N+b2W0AFgA+CeC1betVZyqEEGKU7GqcqbvfBeAu+t8bw+sfAfAjfdZ57jvToZJv63qDFMNZmpUbL3HxlW7Acrk6q7dZZi13gEps9XDPDs3UzRy76Xp6SFrspo3rOZiVK4q+Bd4+S3VRPl5Nad7okl6y5Bqziln+S0p8sZM0CRc4pOmrdnzJTslGHt29c3IlZ2XWMtISbPThZTmq7NjNQhvKIIbmbOJ1e/Oxi+cLZ/M6n0uInyUHkkQ3L0m5y+ZrkqcnhZTLQS/hO6LFods1m/dMBDMkuJmyeRkzu9XMPmpmD5jZDzfM861mdr+Z3Wdmv7qP/RJCCDFedmRA2gk7vzPtkjZhZjdhfUv9De7+KTN7+q73SwghxHhZJyCdnXGz+5B5u6RNfBeAO939UwDg7g/vYb8qcsmEBmD3GBydScK+jJJnvs4oV85aZNZyB5KSbD2ctVkQw3xBEltHKbcqJZfB8nkIo+Dtl2XXKMSCPuhVIuVOw4fHEl+Ubg/oTXK+a5R9OeM3yr4HLaEN8yDtclBE3J8Zu4uDzFqFNJDMWpRgow+6lJZzYStKsHy3sEzKrMV9zbKJ1+1BEq6OVXMJtoWTXBw/nyqspPn8iPMukmsbKKXdPqXT6nmHdTCpPLzovj+rfT0clMxbcVLaxI00z3MBPNfMPmBmHzwKIa4wszvM7F4zu/dzq0d2tLtCCCHGgJsN+jsNxmJAmgG4CcBLsR5A+34ze767fzrOdBQHdQkAnnFwy3bKYAghhBgdjtp0OGb20Zm2pk1gfbf6IXefA/gTM/tjrDvXe7a9MyxnbGWd1QDs7nmZWWgDS6BWzEv7MGmWYDMyV3Dm2M1kXaCUdvl95GXnjl+3JomF9iklTET5bUJtvO/xsQwHQ8QM2SVvI5Zno2c7lXQZs3lZYgzhzRwucIWky1mYNzp7ed+nU3J0h9zpxcDwi/WyUS7O11OUnePQhiJ/l/N2Z53agFIGZik3hmHwlzLL8DF0o5LoQ1sfhz2HNhx0dNf2CojJ5OEdfNftm9MyEw1hHzLv42kTZnaIddrEZZrnN7G+K4WZ3YC17PvgHvZNCCGE2Jid35l2TJu4G8DfN7P7sbaZ/JC7P7brfRNCCDFWzpYBaS/PTDukTTiANxz99YaljlUeGdq4XN0e3XelpBWlmHoAdrJOkn4WRbhB2VZJp8G9WkulYTD9pFmeZdrCFyJRHs1kXYBKlyUZw5mDmY+HkTxZrJfP5ELi4pAGChsIq81yhOck1x+EgAeWBld0QJbeLCNm5dmmibt3TqENsxj+QDLzxMMBSo9VCctsUdptUyPjsryeKN9eTUIbrlSybrNjl0MbouybZfECZYhD9Vginh9paEPRVJ1LWXhLH6fvroJoxsYOs3l3wlgMSEIIIcQxVtfTHTPqTIUQQowO3ZmeYXYlkVjhLO1eqqwq1xazYNOAB3ZZ9sjmHZipmzl268zhjhnD1c6RNBbLo7FUGc5sY2l9QS7pIA9y4MosaM1Oub3R5cnl2eYs8WX5u0kJthmHDYTpWVKObEZveh5L9LWENsTP/YC0ds7jzYjPu+aUzRvbsmCGZRXSwKXugvOXJPF5UYKtbGN3b3TzLimYIWY5c2hDeU3k128Z3lK2WSrdbqdD2db32/5cwobVGSq5rc5UCCHEKDmtAIYhqDMVQggxOiTzjoyNyqx1LHXEEg07dqfJvBmZbNRLHu3BMpFg+wQxLBfdpFx27PYjyOfUUjiGrTnQAShzjimXgSS+snEWypzVLmByi4Z5+5RnYyk3SsQzljWDtDshF6zFgAV2uy85mPX4JQ9NGCrz8h1GEbbAZdaCXJ25d4FS9ubAi2Xi5q1K5sXHG5Wbt1nKza6XLLyl/o6wE+c7ie4l2JrXc55cv2Ph3HemQgghziaKExRCCCE2wBXacPpECYMDHLIghqHSB5dVY3klOniniRTEJZtq2Si4V2nebX2QXR27mazL7W1hFBE+BhGbcmjD8euDHqENZXk2YBEsvMahCWGby2WzxDcnWXeaZL+ydJuVZ+O82Xn4pA/ofURpl7OKJ1FKZaWWrpFVOLCzVXlR9DGExLuKLLSBgxnmkyjz0rHqUWYtSruVtE6fV5RyWaKPDl5+DJCd5/y9ULj6q++MOB8a2y4aujMVQgghNsBNBqRzQ/Ywv5wuf67yL8uM+OuVf9ny9CQbuxmjBjf4NdfVZJTdiQJ5NRxP3nMKj90MNxBzmvWgGGdKBqRJ8x0mF10v4uKScaZZRRneZp+KMnO6jYx3tVe4akwwGbFRZxLu6LgyTnWnWhiQmuME2YyUSXJpUe+k7SrdiV6l2+hlMUa3uWpMmwEptnNkYHybXHEoKkTLOR3zzIi4gfEuMyBtY51jYqk7UyGEEGI462em6kyFEEKIjXDdmZ5/okzTXn0mvE4ixhbVuNLmcaaZdMnmkwyWcrPtZ2NHM5MRm4oyCbgfoRg2x8MVU3nVmPi+FiQBT4O0y1Fy02Bgm9I6D0gSjpVi5vT5HAaZ8yrJvDMamFzECbLhJsYi0qU9Kcbalm2HlA/nxbjbcvtlQfDcaVlUjWEJNhYHz0xFPcaZZpVhagMSy/lxPDF9dqGtrrLU/MjiYJ4YESuTIkJbuZ4+EmzXMahnhbPk5j07eyqEEEKMFN2ZCiGEGB0ODY0ZFXnc1vB5M3qNL8scuiyPxsofWYWKHnpDn8LdfcaOxn2fL5oviKwYeRtlcXAq3h6LPZN7dUL7M0mq8UT5z6fs2A1S9iyPE4xVZQ7oZIpjHqcsYybuXh6vGguJT8HjTMOlvkpK7ACYhco1mZu3iChEPgaVx4uuClcyvcdJrATT7NAFuDJMs2OXx5WyXLtMxpnGNq44lEV88pjy7PvEkjjB+nupe/TgULKx+vvDdtaZmtmtAN6O9Sjrd7j7Wxrm+xYA7wXwd9393myd574zFUIIcTbZRWdqZlMAdwJ4OYCHANxjZpfd/X6a7ykAvg/Ah7qsV89MhRBCjA4HsDQb9NfCiwA84O4PuvtVAO8G8KoT5vsJAD8J4PNd9ld3pgOZZg7dZfO8XEnCkjYOArBEEo7uXpabJslPpj6FuzMX7mLO8le2ryG+j/aHJfJi+9Nml3JdHDzMS4Pp2XkbpbtpFXsX5Gpy6M4Sp2+sKMPthxN2tobIQmeZt7mKDEugkyCBTigVoHDhUnxfLfvGdZJcHBMeWr6zopRbFeMObVmEYubeBUppt5LEw7HjIRZ8XMsAjubgjqothDhkxcCB+vqOTLcUxNB12f0V+N6MDe5MbzCzKMtecvdLR69vBPDx0PYQgBfHhc3shQCe5e6/bWY/1GWD6kyFEEKMDodVQ6p68Ki73zJkQVsHdL8NwGv7LKfOVAghxCjZUWjDJwA8K0w/8+h/13gKgOcB+F1bK1XPAHDZzG7LTEidO1Mze7a7/2mfPT5r1M644J5N5BPr4ajjAdnTZCB37ZA9fl05VBMJaRXeBzt0q3k7OnbbKtyUARPkckyKJGdUDubwo5U/nkL2Jacvu4sLNy+7POPxICl3EUIbZnSwFiz7hs+LwwWm4WCxm5dlzlmQMq9wJnQIZuDQhuI7idVykn1jUMOEHLtx/7gti37j9xG/JK9SxnCcd07vg4unx5CLSrpN3LyLajp+zs3nQPboIysGDvCjoe7fGZm7N/teOhehDbvpTO8BcJOZPQfrTvR2AK+51ujunwFww7VpM/tdAD/Y5ubtcw/9L/kfZvZ1PZYXQgghThV3XwB4PYC7AXwEwHvc/T4ze7OZ3TZ0va13pmb2rQBeCOApZvZVAD7q/ngdhUsAvnroxoUQQoiT2GVog7vfBeAu+t8bG+Z9aZd1dpF5PwDgegDfifVD2b9tZp8G8OcA/qbLRvZN1wHH7Zm6zdpY2sZyT0cpk+dbJCXZWFKaB1lzynpDD0mpa+k0lkozx26WQ9pH7qqJx53kt2TfKG62eC/G2bwxX5aulpjhOiNZd06F5w9nxztRlWeLzlaS7+fkUJ3GEmhGIRLBabuoXLihHBktV8m+FpuozFqcueWj86I4eHfHbpR22b3L2cXRwVtl8yYOXXZfF1JuVYYvOIaTxxtZWAtQnve2pYCY8865SkBy908A+GUz+4/u/gEAMLMvBvBsAP9+t7snhBDiIuKw81nP9FpHevT6MQCP7WSPhBBCCKgE294ZWqJoW9tgh24t4UT3anfJk2WiIjeXtlGKcd1LsGUl0NixOy9KsJXz9nHsxvfJx6pYBzsnSQIttpHI+byJOU1bCF9gt3N0984X5MKd+onzAbVUGCXhKYU2TIM8OaUDy+7eOD1lV3CQRFkCLsMXyL1LH0Lc9wltP66HAx0ySa4KWwjy+YLyduO8WUgDUDp4OfAiTmchDTzNEnCWvxuvkdmCr4FisvgeqB599MrmTdoGhjGMVVo+VzKvEEIIsW8ceb3lsaHOVAghxCjRnek+8G7SxCYScNcB0W3biDGpWTklWzbLmACwiBmuJJ3GkmN9qppVZd4SSatrEANQSl4s5bJ8G8mcz3loQ3nsZpnkSHJxzBWekJs3uns5t3cZyq5VA/15ehqdpVzKrTmnlqXtZSJrcsBBQZ/vpDAvS7lT7/4IIUq5/KUYS6mxlBuDGa4k7l2gdPCylDsP5wuHNNQyfDiu1eON0LbgtuPp6/jarh7bhDZ+NJSd9yOVYEXJ2e1MhRBCnFscJgOSEEIIsSkbBN3vnQvdmfZxvtUOu2a5KysjVpVgi4O+NyjBFh280x6/5rL83T5BDOxkjNIuv6/Cwdwj17gmC20IrlO6HjOXdOXWjJmtLA0G+W9K5eGq6XC+cAm46O6tHLqJu5edvoXT1lgObQ5bYAk2SuQs8y448KEjVQm0GNqQBDMs2ek70LHLDt15IstzwEN5TTSfO9V5XpVibHa4Z8v1bW+ab2gwBM+bBeFsG5bjx8yF7kyFEEKMEwfOZ2iDEEIIsT8MrjvTs8lQWaRNAp4Wg7W57fg1S67sBozzsiM1r7F1TBbSsG4Pa0wcu5msC5SSVxXaEB3MG2TzFtvg43EQXnNKAxGbF+Tmje5ezu2dhbCHVjdvlBGnLEcGl3aS2wuU7t456W3xGBi7bmOZM6fnG8lHwDIvT2esCod5IvOyPJu01e7e4PxlKTc6n5OQBoBKD1al9qIEXDQV5y+f53VJttCWlWcbOKpgPX12Op82dhl0vwvUmQohhBgfrtAGIYQQYiN0ZzpytjUAupQquQRb8zaznFqWaDjEIbpSK0dqkN9Yyo3L1S7gkq6O3bYghsyxW26jeV/qkIZmiZGNCrOo3R6U8zoHXsSSfYmTk3N7o/x3UJXEa87xraT+cLCy3F6gdPdWublRAuXPOc5qdNmT7LsKYQss61oPmbfYBH0+hQQLlnKP9+8qSdkcVBHbKyk3lk5rKcG2CNN8jRQS8Jw+jyRTt8rmzcJKesizCnEYJxeuMxVCCHE2kAFJCCGE2AiTzHvaxDCGbIBxLa2QrDpwkHO1nVTeOX5duWerAIHQRvuT7V5cjvckC5gYGsQAANMgh/WTuxqbqs+r/Gwpwza801mLmzfKx5y9Gt29WW7vnBy6sxnn7ybO3465vQAwD2EMtDuFBMvu5mWWqUuy7yrU96uyeXuENkTzSJbNW7l5s/zhFU+HzyCZd8Euac7fDdMs0cdzgt28xaOPlmuiuH76lFVLrtE+Tt+zhkOhDUIIIcTGyM0rhBBCbIiC7kfEJhmUXedlyYYHZEcJsna2njwfkMtEnDfbr8ZW8zaiBNsniGFauRxPXieQH+c+IQ6T5LjGY8lO3wlJqdNQgo1URCyCjMdOWw9S7oK2vyC5eBoCHthJ2jW3l/d9RnW7YsDDgmXdMOthS3Z4lGQnxnJxnwJ/x3DptPiZVFJulreb5O9WgRdJNm/m2mapP8vmjdfIwZzP8+bQk40eGw2Ufftt4/j1PrN4I5J5hRBCiE1xk8w7Zob+cutjjOFxhFaYjEDznvwaAJwcJtHjUt8gDBv/VxUwHljtJTMZ1YWQu40zbaWI6CubStNReWzS40pJe8sYGUh3tLEw9HRm1NZsSKoKfidRg3MeaxzuFM163DLE1fI5SEXPD8JJWo8rDcaylnMu3n3yEId498nmpBiTyONMryaRgZU5KTF9ZYakujpQeB+sOKTneTldmoya1ZqsbRP6VMkaA+s709Pei+7spVicmd1qZh81swfM7IeT+b7FzNzMbtnHfgkhhBDbYOd3prb++XwngJcDeAjAPWZ22d3vp/meAuD7AHxo1/skhBBi/Ci0oeRFAB5w9wcBwMzeDeBVAO6n+X4CwE8C+KE97BOAzYrw9pFeyqoxJCFFE01VOJy2Gc0nZHaIqpUn6h+biqp97VjthbefmYzq6jNxueZ9qceV9hjj2GOcabEcSbDRVHK1koeD+YXGlXL1mWhIMooBzKIG2QA0i2Niad5YuJsjLhGLbPcxIJGRKY4zbbsEsnGm8UuSTUaLQrq1xjagNB1llWFY1uXi4FGWZ/OYJ9dvvA6yKktAea7XEnB3c9I2xr/3qUxzWuzSgGRmtwJ4O9YXxjvc/S3U/k8BfA/Wp/lfA7iDbwCZfci8NwL4eJh+6Oh/j2NmLwTwLHf/7WxFZnaHmd1rZvd+zh/Z/p4KIYQYDaujFKS+fxlBLX0FgJsBvNrMbqbZftXdn+/ufwfATwF4W9u+nroBydY/098G4LVt87r7JQCXAOAZs1vO0KNpIYQQfXDsLLShVS11978K8z8JHRye++hMPwHgWWH6mUf/u8ZTADwPwO/a2lX4DACXzew2d793D/v3OEMLgLdJJmlU2Kq5jYyd5EJtLg6+TD53dugymWM3SlptY0czx25eUaZ536p4xeL4dK8oUxdxPn59wBLfNM5H0nqQUufs8kzcvX2iBuckCc+mwWnLVWPC/tRFvOObpjGf5OaNUm41ljWs1kiCzp5v1Y7d433ngflXEul2QdNXg427km6XzTJvNSY1RgaSEzt+toccsZmd58n3SZ95h45/H6N02wu3XT0zPUktfTHPZGbfA+ANAA4B/L22le5D5r0HwE1m9hwzOwRwO4DL1xrd/TPufoO7P9vdnw3ggwD23pEKIYQYF6uVDfoDcMO1R4JHf3f03ba73+nu/zmA/xHAj7XNv/M7U3dfmNnrAdyN9U/id7r7fWb2ZgD3uvvlfA1CCCEuGhvKvI+6e9MQyza1lHk3gJ9r2+Benpm6+10A7qL/vbFh3pf2Xf+uBiOXEmT3R7SVJLxqduyWRcZpPYkcWUsKzdvI6FO4u2sQAwAcBAdtLym3R5HkKGdXAROpBNy8DQ50iE7gFbVFdy/H/mXuXnaLRrl2Os3dvFcXHYMa6Mo+LKY4hrBc50E48fjLrKga03JJxGVZyo0uTXbsRmmXQxuu0PuPRb/r4uDBTUz7yp9B3FeuDBOLhVcRm4kzvQ5E6XZubzLKQHTicbUU6070dgCviTOY2U3u/h+OJv9rAP8BLZy6AUkIIYSo8N0Mjemolr7ezF4GYA7gUwD+cdt61ZkKIYQYJbsKbWhTS939+/qu89x3piwBT9Ji4d3bykHWlMWbyJgsE8V562zecjq6eyuHavMmU/JAhe5BDAcUjJDJ17uonsGVLaIEzAETdVYvGueNQRF1qEa33F4gLyQ+KYIYKKSAHLtxkiXgqEJzpm9kRucOO3ZX4YBwAfI+xcEjLBcvYyUWaovSLld7qdy9of0qVS6Ky3JWMk9H2XdJ65nFQvdVXnU8l1nWLSbLayvN5kVjG3PW8nb74DBVjRFCCCE25SwF3aszFUIIMTrcSwVj7JzLzrRrYdttueaywtQAu/go63TZLBWylFq6BWng/8Di4EMdu5msu543biORtLZ2sXQPsahd0s15xLHo+Ywao7s3y+0FyuzkBTl9J4vMzVvuT5SEOdChDFHoXp6NQxtiCTbn4uDhuHJoAxOfd7FcF0McOIhhnoQ29AtmaA7D4GLuMaiB3bzxY+fHANMim7dcLj3vTyEb9yy6gBV0L4QQQmyInpkKIYQQG+CQzHsu6SPLsHQap1n6WRYZvywVNrtO6+CBbk/q22TVro7dtiCGzLEb9yErCdcW4LCMTtv0eLRdkM3zrsJn4FXZu+PpGcm6c5KEpwfHrzkwYBpk1vmk2ekLADaJTlsKdCikXT6wx238a/+A9jUGV7CUOzS0oWqLUm7l9D1uu0KuWw6tKEqnkcx7dT49cT4A1yLnQntzqMYT0jJrx68PrjbLuuvpbi72zOnL023ziv2hzlQIIcT4cA2NEUIIITbCAXgyZn9sqDMN9BkAnTvzWHpJyrVFyXPVw3XKDtUkwKCcr3mdQEvWaMcgBqB0PfaRcvvJVGFZWm5ZHIM8mzfK55W0Ht7Hko7rNAltWNDxiUEApEYiZgQYh4zQ/sT9u0KO3evCQViSXFzm8ZZt7Jichm8wDmkYqiKylLsqsnCbHbtZ3u563hDMwI7doqwau3fpHAjLzii0YVq42Jsd9lXISZatXZ33278DOw+Sr+5MhRBCiE3w+tn2mFFnKoQQYnRsWIJt71zozrRPaEOf3N66PYYCcKDD8WuWPDmb9yDIlfOD5pNs2uM5Qx3aEF73CGLgwezZ+0qPZY99j/OuWNVMHbz0GQS1NJPWOagiul6XXLotcffalI5HLM9G65ks2N17/PqQRFfO8Y1MwzadHkTNppzNe3xAliTztgU1RNLQhrAaHv5QlGAjyTUPZmjO3+WMX85OjtMzOgfiuZ090qlyt9md3/E7pP2xUfN6zhuuO1MhhBBiOA5l8wohhBCb4abQhjGzLYlkqOybST9LLuVG7tEYTHBAUmUtc3ZjqGM3k3V53trBHJfbzsXC7uaSTbYRQhvoGMdjsCLp0NjBkgAAGNZJREFU1umzuzJNpNwYksAZv5SbuygCFfp86M1uXpZgYxhEFdowGXarwF+KUQKuyqyFc6IKYkhKqXHYQszfrdy8JB97aJ/yuV047psfb/DjjOoRypZKD26D097+eeTCdaZCCCHGj0NuXiGEEGJjVDVmRLCcsUrecZ+whWwbdTDC8WuWfpZJObJq3iA55lm0zbQNDu/q2M1k3fV2wmt2yC6jNJbuTme4BF253pbQhiDJ8v5EaXdGbt4YjsHHtZq3Y7m2LNABqEMduhKHGByQm3dVlX07nmZZN5Yn42zgbIA9y7xFaAMdu3g3ksm6ADAPB+jqVXL+Buk2c+8CZRjDQdZWlR48+XXb9NA25jRKue0NB1ZKQBJCCCGGI5lXCCGE2BRXCbYLQebM4+lpkIuXqbOV2pIggiVn83Z0rLZl6k4S2blPEEPm2J0mkla5zvw9RXmykvMLN20uic+udjt2LIdG1+eMpNt83mGBDgA7gbNybd2zeVecCR3WU4U29DAQRzW5Cm2IucZcgm0Z5dmWbN4kmGEZ2jL3LlCe61lIQtUWl0uc8fV6eN7mNuZcSbkJDtOdqRBCCLEpZ6lqzMDRiUIIIYS4xoW+M92XnJJJwqV02ramKLmR/BHDH1hiTOTStCRcjyAGHug+VMplGTojl7b7hDiEkIJEauf3uApyqFPbckpSatAAnRy5XQMdjv5z/PKQDtbVY237gAJm4xCDKVkk2bE7nTYfOy7J1hWWcmPmKpdOKyRgartytUzDyIIZ5ombl0upRac6B5IULvbKmR5fDz+XM4Z+D/UpKTlKXEH3QgghxEbIzXuG6XMHlRb87jH2q4wT5NYsoq75DqFPoeGhJqPsTpTn7XP32WvcaXE3Ttsvjl1+NxWP60FiRuKi6x7GHPL28wozTvN2G4N6tBePv+LowZLm5VaUPclmqVjovBpnOjCasRpnGqZ5YH45PpSrxtDdZxhbynefi3Bc2XDE529RGSYp8s3nR5w+uFo0tcQJ8rXVPWow3nH2uWs9i8alXY0zNbNbAbwdwBTAO9z9LdT+BgDfCWAB4BEA3+Huf5atU89MhRBCjA9fPw4Y8pdhZlMAdwJ4BYCbAbzazG6m2T4M4BZ3/2oA7wXwU227q85UCCHE6Lgm8w75a+FFAB5w9wfd/SqAdwN4VbFt9/e5++eOJj8I4JltK71wMm/xUP5wS+vsIeuyxJnF7lXLZoXESWbsso6T1hPhSLwsBjAtmlyNbW3evz4SdQnLkWGdPcaZLg6bzVtOxyPCUilL9HEcKleYKecrf99eQXnwrgu/f0lVhAd5djVzajt+PaW2dJwpxxkOLQ5O50BhMkokYJZ1uTJMlHZ5LGmc95AfS2TVX3jejhJwJt2up09+fdJ017ZN5h09XsZXbpEbAXw8TD8E4MXJ/K8D8DttK71wnakQQojxs2Foww1mdm+YvuTul/quxMy+HcAtAF7SNq86UyGEEOPDAR9e7/hRd7+loe0TAJ4Vpp959L8CM3sZgB8F8BJ3v9K2wXPZmUapgx2YQ9ezrUoOdVWSML6NdDuWHKNzkN8XxwtG+o0zPX7dJwYwc+zup0hys6TGEuwiHWfK4zw3H4MKlONQeQxqdNrWcn1z9GAeGdi8nqWTzM1VY5IC4NsqDl64eWm3YyxgNnYUKKXdK1dZuj1uy9y7QCnfshM7Xj8cPWnJ44yhkYHDH3WIjtwD4CYzew7WnejtAF4TZzCzFwD4BQC3uvvDXVZ6LjtTIYQQZxvHbp6ZuvvCzF4P4G6sh8a8093vM7M3A7jX3S8DeCuAJwP4DTMDgI+5+23ZetWZCiGEGCW7Cm1w97sA3EX/e2N4/bK+67zQnekmhXX7VXlI4gSTyieV5LiloIZiuR6Fu/sEMUyTefscu670qRozQ7Pse3A1d+U2kQU6tM0bJdhZi6V7EeadV7JulH2pAHj4UprOmtuAUubl+MDBoQ1J1Rh2+hbVXriNHLtR2p1R28GV8Lkucpm3cOxWhd6br99p8ghnV47dC4OKgwshhBCbYwPvTIc92d8MdaZCCCHGh+fmyYzTyPg/951pLv8NX0/Z1ixjrqebpdwoDTnnsi6bJbZM8s1g6ZjJHLtRCsscurzsUOdim/SVf5bdq8ZE2ZedvvEzYKdt10AHhnN8I3wOZCzIzXslaGIzyt+dHYQgBmpjc3GUeReJe5djgzNJjiPeytAGXk8IbSA5lkMborQ748owYXpKEjDLvNHBy23Rxc0qfHT3toc2dHOxt2d978INPz4M7UE2Y+Lcd6ZCCCHOIG5V4tSYUWcqhBBilPDY7jGjzjQhK647VKZJM2x58HyyjaGuym1l6rZL293n7drWj7jNNkm8+TPgQfolSbF2IhZs51CAuOysh8xbcyxlLhI37wFn81ahDSfvG1BKwm0fVTZGsAht4LJqQZ5l6Zgdu1HarR26x/NmxcDX8zavJ4alVKENRTZvua/buibE2UCdqRBCiNFhXhckGDPqTIUQQowSGZDOIZu57+JrlnfCAPkkt/foP43byJytfdy0QzN1+zgOd1FuKn//5CSleWMG8gHJePPDbqEanGfLz3p4vSVx2T55u7w/YS3k2J0HN+8VzpCl86wIbaDdGapGZo5dlnLj2MKqdBpLsIlj9yApq5ZJuVUGcyLl5nnVaKTfo6Dm9ZxnzM9WTrE6UyGEEKNkaGjDaaDOVAghxOgwr9WxMbOXztTMbgXwdqwT+t/h7m+h9jcA+E6sgyseAfAd7v5nu9iXftJLsyM0a+uz/UwC5gzZWMbKeaB9x5OOHbrM0EzdfvM2b3+4rNPshG4L6ogOzazsXb7Nbf2C5nOp/KAnq1ASbkXv+aD5w43j9bg83Iql3Gm30IY+8B1GUdaMZedVs+TKGbuZYzdKu+ygTkuwcSm3q81ycfk4o7vDfRMujux7tsaZ8gOarWNmUwB3AngFgJsBvNrMbqbZPgzgFnf/agDvBfBTu94vIYQQI8bXPxyG/J0GO+9MAbwIwAPu/qC7XwXwbgCvijO4+/vc/XNHkx/EuvK5EEIIcSbYh8x7I4CPh+mHALw4mf91AH7npAYzuwPAHQDwlMmXb7xjQ52jbW25s7VZLua2WjoM+arNu5OySaZuPzdv8z7sJlu0qxybk+1P/fwmk/p3Jftm6z3+bczn0iJm8/IjApLBYxACS8IRloezYQyZXJdJuVk5NF5v5fQN0m6bmzeXcsN6aF9nV0+e76R9H3ptXVQMMiANxsy+HcAtAF5yUru7XwJwCQCeMbvlNKrsCCGE2AcyIFV8AsCzwvQzj/5XYGYvA/CjAF7i7lf2sF9CCCFGiqrG1NwD4CYzew7WnejtAF4TZzCzFwD4BQC3uvvDe9injennbG0e+N+1VNmaZimI81Wz/Wnafr0/wx26Q8tNDSUPsWiWyxnOXo3u3iqnN7S1lc/bjezb/PksDmjrQS6bsZt3mrl7y20sY6BDy55mUXBlEAIHh0QXbu6QnSWZutG5ft3nKdOXHbtREk5KD2Zl1aLke9K+DpVys4zwat4zdCfXikIbStx9YWavB3A31tffO939PjN7M4B73f0ygLcCeDKA3zAzAPiYu9+2630TQggxXtqG8o2JvTwzdfe7ANxF/3tjeP2yfeyHEEKIs8E6tEF3pmeSzSTYruvJJODucmS9zW4nXZsknc/bfT1dj93pyFLDZF8OdCilwhZZ97BZnh0KH7so7bILMkq7qyllFZPsuwz6rVPbrMe+Z07MIhxk1Xzu8Do4fKGQWRPHbibrcjvPG/c1k3I3CW0Yeh2cK1n3BM7S+1NnKoQQYnSY50OqxsY+QhuEEEKIc43uTAfSRwLu6oKt19MsG2Xu3Xr725dyN3vP+/i12b1cXZQ9WcqN9HL6ZrJvtY3uxyP+Uo9BDEBp1pgcsEP3eDkuhzahn9RRAmVJeChV+bqOEnDm0AXKHF2+iymk20TW7TPvZg735rZsub7t5wkO9Bgz6kyFEEKMDzcZkIQQQohNMD9bd+HqTAfS1aHL7Vnb7CrlqR5me9Bduu3axu27KqvW9QJpmy8vrTYsq5flv7iNzM2bO315m3Q8pt2dvtl7LjNkKWwhbGNGsm4le8d2kjyHkub2VqENx6/5fVSyb/Keo1ybybq8zex8ZTdvKQGXbfnjle7XxFkKLtg2u0pA6lAW9BsB/HMAXw3gdnd/b9s61ZkKIYQYHw7YDn5IhLKgL8e68Mo9ZnbZ3e8Ps30MwGsB/GDX9aozTSjuxDjSK7lr7GfqiVP865lj35rX05XTGGfatt6ubX3Ii4MPG8/bx4DUZ0zqctXcxndbcb18zDMjVXZc2WS07GFuG8o0MbNlxcH5eBR341wcvLijzMeZZrGA+d3nya+7THddzy7ga2KMcqphZ0H3j5cFBQAzu1YW9PHO1N3/9Kit872xOlMhhBDjY7NnpjeY2b1h+tJR1TGgf1nQTqgzFUIIMToMGz0vftTdb9ni7rRyoTvTPtJttezg8WW5WalkeLxguY3dSLld19M2b9e2bN7cjNRGFgvYTB8DUmoc6mgqYpYs1a26ScAs4055XGUytjQrFt4nrSYzDhXrTCq4rJftNm9d1DszGfHjljhfZsrrHid4GrJqn8dEcd7Nrq1R0qksaF8udGcqhBBipPjOqsa0lgUdguIEhRBCjI5rBqQhfxnuvgBwrSzoRwC851pZUDO7DQDM7O+a2UMA/iGAXzCz+9r29+zemVqz/DDU9dpvrFdzXF2fbfRb9ngfuo437NPG7UMrwfSdt1xu6DOS4Q7UUi5ulmc5wjGOA65l3aHz5uspi3yD5j35NVBKwrNsXGkLfeIF0/MwueOIwyH4izF79JBGDfaSeZu3kVeN6b6vqQTcoxj4uWaHxcE7lAW9B2v5tzNntzMVQghxflECkhBCCLEZazfvae9Fd9SZJgx3pHZ37Jbz9qkE03nWrUX9bStOsOty/egeH9hWRabrNmKoBq9jcdh93lI65aAOjh6Mbc3hD9Xnk0rJaIQdw+z8zcieXXU9X1i6rWXfroEK+XmeOXbLNl6u2/ZPmu7KaYQ4pPPuq9fYocy7C9SZCiGEGB1n7c5Ubl4hhBBiQ3Rn2pE+ku/wAINciupaEHwTWXVslWG2Q5tUNNQJ3CzPsuszyrW1zNtdgo3zZhIwE9fTJt16+Im9rS+IbLxgJt0yQ6XctnO5q1ybuYB3FXLStuy5RQYkIYQQYkPUmQohhBCbYTAZkPYNu8v28WsmDqxmeS2Xe5qdpW0u0yy0Yais2k++bpbC2pbt2rYL2t27ze+rlFWb15OFPbSvp9v2ed58f1i67hb40aV9CNtzxjcvu8l6uufvNm8/k4A32dcLi+5MhRBCiM0wdaZCCCHE5qgzPSdsa5B15tjNBvP32cbQeYc6EMcm8w4PYiiXzSTgzfJ/u7l580CHXMotl+0ezLCPElu7KsnXPbShebm2bXaVeTdy9W8pj7fIJU/W2faYaAxl10yhDUIIIcTmnKU7U4U2CCGEEBty7u9MN5FOs1JupbyUyW1tv666OX1PWm8Tu8ji7Ttv17azBB+77pJrW8m8roEO3aXcXNYdLlcPpY9ct71zsntubtdtbBIM0ecxSbaermQZun2+W05N8pUBSQghhNgMuXmFEEKILaDOdA+4na7jrJBsyDWXZaTm5FLYLko47UPK3ZZT8TTInb5xvuGhDflyzedEmyTcvNzpOySHhhZs6/zcxaOQXjIvf2f0CES5KMjNK4QQQmyBs/TDQp2pEEKI8aFnpqdP5sIdSh/322mfAJsNgu+4nhbp9vTlmWbJk/c9uh6Hfs5tjxwyx265vdwVnLd1z989bYZmSfdp34dcnLWd9vdAxtBs75OW3RUyIAkhhBBb4Cx1pgptEEIIITbkXNyZVpLFltyjmYyXDeQemuG6Kyl5Wy7cbQ283w99ZOZun1efUI+hg+DbnLb7CGPYxWMSpuu5tKtwg9O4fsrlTvsxSMkYHwtI5hVCCCE2xc/WsDp1pkIIIUbJ2O7gM9SZdqSPxFdLjNuRfbex3CYu3PObv9vtPbOMGo8l56AOlYTbjmMux3X74tlE0tuVO36sy1Xr2eBOaVvlDYdSuNZ7vI8s43eX7FLmNbNbAbwdwBTAO9z9LdR+HYBfBvC1AB4D8G3u/qfZOmVAEkIIMUomy2F/GWY2BXAngFcAuBnAq83sZprtdQA+5e5fAeBnAPxk674OeYNCCCHELrl2Z7rtzhTAiwA84O4PuvtVAO8G8Cqa51UAfuno9XsBfLOZpdLPmZZ5x5LNy9TO3+Yya1vbn4HOwYteVg3oHr6Q5e/W0m2zJMxsS0bbx/UwBrm0cZ0jf76272smG4FQPZYYo9FnM5n3BjO7N0xfcvdLR69vBPDx0PYQgBfT8o/P4+4LM/sMgC8G8GjTBs90ZyqEEEKcwKPufss+N3guOtOxjZHKf001V6RoGye4jbF5m8Wzddp8K31i+MZMPg65+13S6oze5Z+2Gek02Nb5uqsxsUPHCPdRR/Z5ze7o3PgEgGeF6Wce/e+keR4ysxmAp2JtRGpEz0yFEEKMjh0+M70HwE1m9hwzOwRwO4DLNM9lAP/46PV/A+Dfunt6t3Mu7kyFEEKcM3YU2nD0DPT1AO7GemjMO939PjN7M4B73f0ygP8DwK+Y2QMAPol1h5uyl850F2N6MjLJYpMxmEPpKr1sYqDYRxWMbe1PNt8uZLQ+6xy6P5tU3diHzLmLbYxd5h16Lg2N8Rz7I4vBEvApvS/D7s4Nd78LwF30vzeG158H8A/7rHPnMu+uxvQIIYQ4x+xO5t0J+3hmupMxPUIIIc43Z6kz3YfMu7UxPWZ2B4A7jiav/OxD9kc72eOzzw1IxkNdcHRsmtGxaUbHppm/vYuV/gV+/+43wW4YuPjeP6szZUA6GnR7CQDM7N59jyM6K+jYNKNj04yOTTM6Ns1QOMLWcPdbd7HeXbEPmbfPmB50HdMjhBBCjIV9dKY7GdMjhBBCjIWdy7y7GtODI7lXnIiOTTM6Ns3o2DSjY9OMjg0A0w2gEEIIsRmKExRCCCE2RJ2pEEIIsSGj70zN7FYz+6iZPWBmP3xC+3Vm9utH7R8ys2fvfy9Phw7H5g1mdr+Z/aGZ/Rsz+1unsZ+nQduxCfN9i5m5mV2YYQ9djo2ZfevRuXOfmf3qvvfxtOhwTX25mb3PzD58dF298jT2c9+Y2TvN7GGzk8f225qfPTpuf2hmL9z3Pp467j7aP6wNS/8RwH8G4BDAHwC4meb5bgA/f/T6dgC/ftr7PaJj800Annj0+p/p2FTzPQXA+wF8EMAtp73fYzk2AG4C8GEAX3g0/fTT3u8RHZtLAP7Z0eubAfzpae/3no7NNwJ4IYA/amh/JYDfwTpS9+sAfOi093nff2O/M1UUYTOtx8bd3+funzua/CDWY3wvAl3OGwD4CaxzoD+/z507Zbocm+8CcKe7fwoA3P3hPe/jadHl2DiALzh6/VQAf77H/Ts13P39WI+0aOJVAH7Z13wQwNPM7Mv2s3fjYOyd6UlRhDc2zePuCwDXogjPO12OTeR1WP9yvAi0HpsjGepZ7v7b+9yxEdDlvHkugOea2QfM7INHVZ8uAl2OzZsAfLuZPYR11ZHv3c+ujZ6+30fnjjMVJyiGYWbfDuAWAC857X0ZA2Y2AfA2AK895V0ZKzOspd6XYq1mvN/Mnu/unz7VvRoHrwbwLnf/aTP7eqzHxz/P3VenvWPidBn7namiCJvpcmxgZi8D8KMAbnP3K3vat9Om7dg8BcDzAPyumf0p1s94Ll8QE1KX8+YhAJfdfe7ufwLgj7HuXM87XY7N6wC8BwDc/fcAXI91CP5Fp9P30Xlm7J2pogibaT02ZvYCAL+AdUd6UZ57AS3Hxt0/4+43uPuz3f3ZWD9Pvs3ddxLYPTK6XFO/ifVdKczsBqxl3wf3uZOnRJdj8zEA3wwAZvZVWHemj+x1L8fJZQD/6MjV+3UAPuPuf3HaO7VPRi3z+u6iCM88HY/NWwE8GcBvHHmyPubut53aTu+JjsfmQtLx2NwN4O+b2f0AlgB+yN3PvdrT8dj8AIBfNLPvx9qM9NqL8OPdzH4N6x9YNxw9L/5xAAcA4O4/j/Xz41cCeADA5wD8k9PZ09NDcYJCCCHEhoxd5hVCCCFGjzpTIYQQYkPUmQohhBAbos5UCCGE2BB1pkIIIcSGqDMVQgghNkSdqRBCCLEh6kyF2DNH9TBffvT6fzKzf3Ha+ySE2IxRJyAJcU75cQBvNrOnA3gBgHOfSiXEeUcJSEKcAmb277COenypu3/2tPdHCLEZknmF2DNm9nwAXwbgqjpSIc4H6kyF2CNm9mUA/i8ArwLw1xeo8LYQ5xp1pkLsCTN7IoB/CeAH3P0jAH4C6+enQogzjp6ZCiGEEBuiO1MhhBBiQ9SZCiGEEBuizlQIIYTYEHWmQgghxIaoMxVCCCE2RJ2pEEIIsSHqTIUQQogN+f8BfFFaZDhbcZwAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#Compute the Error\n",
        "error_s = np.linalg.norm(S_test - S_pred) / np.linalg.norm(S_test)\n",
        "print(error_s)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WkDnT1Mj5YWm",
        "outputId": "15e9fe11-5896-4ae7-d45c-23b71ea240bd"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0.006711102955147286\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# References\n",
        "\n",
        "[1] Lu, L., Jin, P., & Karniadakis, G. E. (2019). Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint arXiv:1910.03193.\n",
        "\n",
        "[2] Wang, S., Wang, H., & Perdikaris, P. (2021). Learning the solution operator of parametric partial differential equations with physics-informed DeepONets. Science advances, 7(40), eabi8605."
      ],
      "metadata": {
        "id": "IZgbjRR155tS"
      }
    }
  ]
}
